{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 14\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from previous chapters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(beta, gamma):\n",
    "    \"\"\"Make a system object for the SIR model.\n",
    "    \n",
    "    beta: contact rate in days\n",
    "    gamma: recovery rate in days\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(S=89, I=1, R=0)\n",
    "    init /= np.sum(init)\n",
    "\n",
    "    t0 = 0\n",
    "    t_end = 7 * 14\n",
    "\n",
    "    return System(init=init, t0=t0, t_end=t_end,\n",
    "                  beta=beta, gamma=gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the SIR model.\n",
    "    \n",
    "    state: State (s, i, r)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (sir)\n",
    "    \"\"\"\n",
    "    s, i, r = state\n",
    "\n",
    "    infected = system.beta * i * s    \n",
    "    recovered = system.gamma * i\n",
    "    \n",
    "    s -= infected\n",
    "    i += infected - recovered\n",
    "    r += recovered\n",
    "    \n",
    "    return State(S=s, I=i, R=r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    init, t0, t_end = system.init, system.t0, system.t_end\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t0] = init\n",
    "    \n",
    "    for t in linrange(t0, t_end):\n",
    "        frame.row[t+1] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_total_infected(results):\n",
    "    \"\"\"Fraction of population infected during the simulation.\n",
    "    \n",
    "    results: DataFrame with columns S, I, R\n",
    "    \n",
    "    returns: fraction of population\n",
    "    \"\"\"\n",
    "    return get_first_value(results.S) - get_last_value(results.S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_beta(beta_array, gamma):\n",
    "    \"\"\"Sweep a range of values for beta.\n",
    "    \n",
    "    beta_array: array of beta values\n",
    "    gamma: recovery rate\n",
    "    \n",
    "    returns: SweepSeries that maps from beta to total infected\n",
    "    \"\"\"\n",
    "    sweep = SweepSeries()\n",
    "    for beta in beta_array:\n",
    "        system = make_system(beta, gamma)\n",
    "        results = run_simulation(system, update_func)\n",
    "        sweep[system.beta] = calc_total_infected(results)\n",
    "    return sweep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_parameters(beta_array, gamma_array):\n",
    "    \"\"\"Sweep a range of values for beta and gamma.\n",
    "    \n",
    "    beta_array: array of infection rates\n",
    "    gamma_array: array of recovery rates\n",
    "    \n",
    "    returns: SweepFrame with one row for each beta\n",
    "             and one column for each gamma\n",
    "    \"\"\"\n",
    "    frame = SweepFrame(columns=gamma_array)\n",
    "    for gamma in gamma_array:\n",
    "        frame[gamma] = sweep_beta(beta_array, gamma)\n",
    "    return frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contact number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the `SweepFrame` from the previous chapter, with one row for each value of `beta` and one column for each value of `gamma`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0.2</th>\n",
       "      <th>0.4</th>\n",
       "      <th>0.6</th>\n",
       "      <th>0.8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.1</th>\n",
       "      <td>0.010756</td>\n",
       "      <td>0.003642</td>\n",
       "      <td>0.002191</td>\n",
       "      <td>0.001567</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.2</th>\n",
       "      <td>0.118984</td>\n",
       "      <td>0.010763</td>\n",
       "      <td>0.005447</td>\n",
       "      <td>0.003644</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.3</th>\n",
       "      <td>0.589095</td>\n",
       "      <td>0.030185</td>\n",
       "      <td>0.010771</td>\n",
       "      <td>0.006526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.4</th>\n",
       "      <td>0.801339</td>\n",
       "      <td>0.131563</td>\n",
       "      <td>0.020917</td>\n",
       "      <td>0.010780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.5</th>\n",
       "      <td>0.896577</td>\n",
       "      <td>0.396409</td>\n",
       "      <td>0.046140</td>\n",
       "      <td>0.017640</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0.2       0.4       0.6       0.8\n",
       "0.1  0.010756  0.003642  0.002191  0.001567\n",
       "0.2  0.118984  0.010763  0.005447  0.003644\n",
       "0.3  0.589095  0.030185  0.010771  0.006526\n",
       "0.4  0.801339  0.131563  0.020917  0.010780\n",
       "0.5  0.896577  0.396409  0.046140  0.017640"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_array = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 , 1.1]\n",
    "gamma_array = [0.2, 0.4, 0.6, 0.8]\n",
    "frame = sweep_parameters(beta_array, gamma_array)\n",
    "frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11, 4)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "frame.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following loop shows how we can loop through the columns and rows of the `SweepFrame`.  With 11 rows and 4 columns, there are 44 elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1 0.2 0.010756340768063644\n",
      "0.2 0.2 0.11898421353185373\n",
      "0.3 0.2 0.5890954199973404\n",
      "0.4 0.2 0.8013385277185551\n",
      "0.5 0.2 0.8965769637207062\n",
      "0.6 0.2 0.942929291399791\n",
      "0.7 0.2 0.966299311298026\n",
      "0.8 0.2 0.9781518959989762\n",
      "0.9 0.2 0.9840568957948106\n",
      "1.0 0.2 0.9868823507202488\n",
      "1.1 0.2 0.988148177093735\n",
      "0.1 0.4 0.0036416926514175607\n",
      "0.2 0.4 0.010763463373360094\n",
      "0.3 0.4 0.030184952469116566\n",
      "0.4 0.4 0.131562924303259\n",
      "0.5 0.4 0.3964094037932606\n",
      "0.6 0.4 0.5979016626615987\n",
      "0.7 0.4 0.7284704154876106\n",
      "0.8 0.4 0.8144604459153759\n",
      "0.9 0.4 0.8722697237137128\n",
      "1.0 0.4 0.9116692168795855\n",
      "1.1 0.4 0.9386802509510287\n",
      "0.1 0.6 0.002190722188881611\n",
      "0.2 0.6 0.005446688837466351\n",
      "0.3 0.6 0.010771139974975585\n",
      "0.4 0.6 0.020916599304195316\n",
      "0.5 0.6 0.04614035896610047\n",
      "0.6 0.6 0.13288938996079536\n",
      "0.7 0.6 0.3118432512847451\n",
      "0.8 0.6 0.47832565854255393\n",
      "0.9 0.6 0.605687582114665\n",
      "1.0 0.6 0.7014254793376209\n",
      "1.1 0.6 0.7738176405451065\n",
      "0.1 0.8 0.0015665254038139675\n",
      "0.2 0.8 0.003643953969662994\n",
      "0.3 0.8 0.006526163529085194\n",
      "0.4 0.8 0.010779807499500693\n",
      "0.5 0.8 0.017639902596349066\n",
      "0.6 0.8 0.030291868201986594\n",
      "0.7 0.8 0.05882382948158804\n",
      "0.8 0.8 0.13358889291095588\n",
      "0.9 0.8 0.2668895539427739\n",
      "1.0 0.8 0.40375121210421994\n",
      "1.1 0.8 0.519583469821867\n"
     ]
    }
   ],
   "source": [
    "for gamma in frame.columns:\n",
    "    column = frame[gamma]\n",
    "    for beta in column.index:\n",
    "        frac_infected = column[beta]\n",
    "        print(beta, gamma, frac_infected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can wrap that loop in a function and plot the results.  For each element of the `SweepFrame`, we have `beta`, `gamma`, and `frac_infected`, and we plot `beta/gamma` on the x-axis and `frac_infected` on the y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sweep_frame(frame):\n",
    "    \"\"\"Plot the values from a SweepFrame.\n",
    "    \n",
    "    For each (beta, gamma), compute the contact number,\n",
    "    beta/gamma\n",
    "    \n",
    "    frame: SweepFrame with one row per beta, one column per gamma\n",
    "    \"\"\"\n",
    "    for gamma in frame.columns:\n",
    "        column = frame[gamma]\n",
    "        for beta in column.index:\n",
    "            frac_infected = column[beta]\n",
    "            plot(beta/gamma, frac_infected, 'ro')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap14-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhcVbnv8W8nJEVCg2AGQSEyXHlDEFBG8XAS0IOIciEa7kVtgshwGMIUBUVQlPF4QCBAgmEUjcHohQsowiEchgQVMUEGIeGFDBhQYieAkrEz9flj7Up2V6oru9NVe++u+n2epx+6Vu3e++2qUG+vtddab1N7ezsiIiJ50yvrAERERMpRghIRkVxSghIRkVxSghIRkVzaIusAqs3MCsABwFvA2ozDERGRynoDOwAz3L0t/kTdJShCcnoq6yBERKRL/hX4bbyhHhPUWwCTJ09m++23zzoWERGpYOHChbS0tED02R1XjwlqLcD222/PjjvumHUsIiKSzEa3ZDRJQkREckkJSkREcimTIT4zOxB40N0Hd/L8EOAO4BNAK3C2uz+UYogiIpKxVHtQZtZkZqcAU4G+FQ6dArwIDABOBaaY2a4phCgiIjmR9hDfpcAZwBWdHWBmuwP7A5e4+yp3fxz4FXByOiGKiEgepD3EN9HdLzGzQyscMwxY4O7LYm2vAAfWNDIRkR6mddp0FkyaTNvitykMHMCQ0S0MHjG8bq6daoJy978lOKwZWF7SthzoX/2IRES6J6sk0TptOnMnTGRdW9h8oW3RYuZOmAhQ8+unde08zuJbBvQraesPLM0gFhHpIVqnTWfmKafxu5HHMvOU02idNj2Va86dMJG2RYuhvX39B3Ua114wafL6BFG0rq2NBZMm182185igZgFDzCyepIZG7SIiG8kqUWSZJNoWv92l9p547dwlKHd34AXgSjMrmNlhwDHA3dlGJiKbkkUvBrJLFFkmicLAAV1q74nXzkWCMrMWM4sP4Y0C9iCsgbodONndX8okOBFJJMvhrqwSRZZJYsjoFnoVCh3aehUKDBndUjfXzmShrrs/CWwbezwZmBx7/AZwZPqRicjmqtSLqfVN+8LAASExlmmvpSGjWzpMFoD0kkTxNc1igkZa167HzWJFJANZDndllSiyTBLF66d1rSyurQQlIlWRVS8Gsu9NZJUk6p0SlEidyWpdTpbDXaBEUY+UoERyaHOTTJaLN7Me7pL6owQlkjPdSTJZTlQA9WKkunIxzVxENujOmp4sJyqIVJsSlEjOdCfJZLkuR6TalKBEcqY7SSbLxZsi1aYEJZIz3Ukyg0cMZ7cxp1MYNBCamigMGshuY07XfSHpkTRJQqQGujPVu7uz4TRRQeqFEpRIlVVjqreSjIiG+ESqLssSDCL1RAlKpMo01VukOpSgRKpsi+bmLrWLSHlKUCJV1k57l9pFpDwlKJEqW7t0WZfaRaQ8JSiRKtNuDiLVoQQl0onWadOZecpp/G7kscw85bTEpcu1m4NIdWgdlEgZ3VnLpLITItWhBCVSRnfLVmihrUj3aYhPpAytZRLJnhKUSBma6CCSPSUokTI00UEke7oHJVKGJjqIZE8JSqQTmuggki0N8YmISC4pQYmISC4pQYmISC4pQYmISC5pkoTUrdZp0zULT6QHU4KSutSdvfREJB9STVBmtg8wEdgbmAec5O4zyhxnwI+AfYElwER3vzLNWKVn6+5eeiKSvdTuQZlZX+AB4BfAtsCVwFQz26bM4ZOB/wbeD3wKOMfMjk4rVun5tJeeSM+X5iSJQ4E+7j7O3Ve7+xTgZeC4Msda9N8moD36WplKlFIXtmhu7lK7iORPmglqGDC7pO0VYK8yx14OXAq0Aa8Bd7v71NqGJ/WknfYutYtI/qSZoJqB5SVty4H+ZY5tB74R/czHgC+a2cm1DU/qydqly7rULiL5k+YkiWVAv5K2/sDSeIOZ7Q+MdfcPRk0vmNnVwJnAHTWPUupCYeAA2hYtLtsuIj1Dmj2oWWy4t1Q0NGqP2wnoa2ZNsbY1wOoaxiZ1RuUyRHq+NHtQTwBNZjYWGA+MIkw3v6/kuN8BvYFLzewyYBfgfGBCirFKjiVZgKtyGSI9X2oJyt1XmdmRhHVQlwGvAyPdfZGZtQC3uHuzu7dGx10DnAO8C9wK3JRWrJJfXVmAq3IZIj1bqgt13f0l4JAy7ZMJa5+Kj/8A/GuKoUkPoQW4Io1Dm8VKj6IFuCKNQwlKepTOZuFpdp5I/VGCkh5Fs/NEGod2M5ceRbPzRBqHEpT0OJqdJ9IYNMQnIiK5pB6U5Jqq4oo0rk4TlJkl/hRw9+nVCUdkA1XFFWlslXpQT5Y8bifUZ1oHrAX6RN+vovyO5CLdokW5Io2t0j2orWNfXwP+DBwMFNy9QCiDMQMYW+sgpTFpUa5IY+u0B+Xu6wvnmNnlwCh3nxF7/kUzGwM8DNxS0yilIalkhkhjSzqLbxvCDuOltkYTLaRGtChXpLElTS73AD+OSmU8R7gXdRBwHfDTGsUmDabcjL3dxpyuWXwiDSppgjob+BHwQOxnVgO3AxfUIC5pMJ3N2NttzOnsf7tGkEUaUaIE5e4rgBPN7Gw2VMV9xd2XVvgxkcQ0Y09ESiXeScLMtgZGR1+vA8PNbJcaxSUNRjP2RKRUogRlZrsDTqhwewZh0sSXgRfM7JO1C08ahcpoiEippD2oG4BfuvtQoA3A3UcDkwil2UW6RTP2RKRU0kkSBwPnlmm/DniheuFIo1IZDREplTRBrQAGA6+WtH8EeK+qEUnDUhkNEYlLOsT3E2CCmR0cPR5kZkcTpp5PrklkIiLS0JL2oL5D2Cz2caAAPA2sAW4GLq5NaCIi0sgS9aDcfY27fxvYDtgL+DiwnbufB2xbw/hERKRBJepBmdlaYHt3XwS8HGsfAswCmmsTnoiINKpKBQu/DHwhetgE3G5mbSWHfRh4p0axiYhIA6s0xPcosBQolt1YEX1f/FoKPAOMrGWAIiLSmCrVg1oMnARgZq8D17j78uLzZraFu6+pdYAiItKYkk4zvx6YaGYXxdrmmdkdZqZy7yIiUnVJE9RNwJ7Af8XaRgN7A9dWOygREZGkCeoo4Gvu/qdig7tPA04DRtUiMBERaWyJy20AW3bS3rcagYiIiMQlTVC/Jmx19NFig5ntQRj6e6gWgYmISGNLutXRWOB+4MVoLVQ7YcujqYRy8ImY2T7ARMK9q3nASe4+o8xxWxOS39HRte4BznL31UmvJfnROm26dikXkS5LWvL9XWCEmQ0DhgGrgFfd/ZWkFzKzvsADwDhgOOHe1VQz+7C7l+6IfifQB9iZMLT4MHABcFXS60k+tE6bzpwbJ9C+JqxIaFu0mDk3TgBQkhKRirpS8r0XIWHsBEwHtjGzbbpwrUOBPu4+zt1Xu/sUwrZJx5VcZwfgGOBUd3/P3Vujx9o1vQeaf9ud65NTUfuaNcy/7c6MIhKRniJpyfcdgOeBKYQKuu8HLgRmR/eikhgGzC5pe4Ww+Wzcx4EFQIuZzTOzN4AxwF8TXkdyZM2SJV1qFxEp6krJ91nAIMKWRwDHA3+InkuiGVhe0rYcKF3o+35CT+2jhHtVIwj3or6Z8DoiIlIHkiaoTwGXufv6zWKjbY8uBg5KeI5lQL+Stv6EPf3i2oDewDfcfam7zyOUlv9iwuuIiEgdSJqgmii/DmogYcJEErMAK2kbGrXHFSdexOtMJZ1tKCIidSJpgroXuDa6F9UOYGZ7AxMIM/OSeAJoMrOxZtbHzL5EGMK7L36Qu/8ZmAlcb2ZbmdmHCdPcf57wOpIjhUEDu9QuIlKUNEF9HVhImKjQTOj1PAfMj57bJHdfBRxJmF7+DmF4cKS7LzKzFjOLD/V9DlhJWCs1kw3T06WHGTK6hV6FQoe2XoUCQ0a3ZBSRiPQUTe3t7WWfMLN/B+5296Wxtl2BPQhDbrPd/dVUouwCM9sZmP/YY4+x4447Zh2OoIW6ItK5N998k09/+tMAu7j76/HnKt3bGUdYILs0VvJ9HqFXI5LY4BHDlZBEpMsqJai/EWpA/ZEwSeKCkmG49dz9sloEJyIijatSgjoRuAT4PGFixOFAuQq67YASlIiIVFWlku+/BT4DYGbzgcOjMvAiIiI1l3Sz2F0AzKxP9DNNJc+X7hAhIiLSLYkSlJkdBNxK2H4orokwxNe7ynGJiEiDS7pDww3AP4GRQGlpDBERkapLmqD2Aj4R7fIgIiJSc0kT1Gzgg4ASlGxEC3FFpBaSJqibgNvM7CbgNUo2iHX3h6odmPQMrdOmM3fCRNa1hY3u2xYtZu6EiYAq5opI9yRNUD+O/vufZZ7TJIkGtmDS5PXJqWhdWxsLJk1WghKRbkk6zTxxaXhpLG2L3+5Su4hIUp0mKDPrX1zfZGalVW870DqoxlUYOIC2RRuv3y4MHJBBNCJSTyr1jJaY2eDo+6XAkjJfxXZpUCqnISK1UmmI71OEuk0Ah6UQi/RAxftMmsUnItVWaS++aeW+FymlchoiUgtJZ/GJdKC1TyJSa0pQ0mVa+yQiadD0cemySmufRESqRQlKukxrn0QkDUnLbQwCvgPsB/Rh43pQB1Y/NMkrrX0SkTQkvQd1B3AQ8DNUbqPhDRnd0uEeFGjtk4hUX9IENRw4RtPNBbT2SUTSkTRBvUsoWCgCaO2TiNRe0gR1CTDezMZSvtyG9uITEZGqSpqgrgO2Bf7QyfMqtyEiIlWVNEEdW9MoRERESiStBzUNwMz6AR8hrJ+a6+7ayVxERGoi6Tqo3sBVwLlsWAe1yszuAs5y9zU1i1BERBpS0p0krgSOB04Adoq+TgA+R5hAISIiUlVJ70GdAJzi7g/F2n5pZkuAW1GSEhGRKkvag2oG5pRpnwcMrF44IiIiQdIe1AxgDOEeVNxZwLNJL2Zm+wATgb0Jye0kd59R4fg+hKntv3b37ye9jlSXaj+JSBaSJqhvAU+a2aFsWAv1CWBn4LNJTmBmfYEHgHGErZNGAVPN7MPu3tn+flcAHwN+nTBOqbLWadOZc+ME2teEeTBtixYz58YJgGo/iUhtJRric/eZwL7Ao4QJEgMJSWOouz+d8FqHAn3cfZy7r3b3KcDLwHHlDo6S4eHAIwnPLzUw/7Y71yenovY1a5h/250ZRSQijSJxRV13fxU4vxvXGgbMLml7Bdir9EAz2w64DfgCoRclGVmzpPxSt87aRUSqpdMEZWZ/BI5w93fNbAbQ3tmxCetBNQOle/YtB/qXOXYicLO7v2RmCU4tIiL1plIP6jdAseDPg1W41jKgX0lbf2BpvMHMTiQMIY6rwjWlm3pv3czaJUvLtouI1FKnCcrdL409fAJ42t1Xx48xswJhsW4Ss4CxJW1DgZ+WtH0ZOBB4N+o9bQV81sz2d/ejEl5LqmTXU0/mtRvGw9q1Gxp792bXU0/OLigRaQhJ70E9AWwPLCpp3xW4m417Rp2doykq2TGeMItvb+C++EHufkT8sZndDzyvaebZUHFCEclKpXtQZwDFXlQTMMvMSu9DNQPPJbmQu68ysyMJ95cuA14HRrr7IjNrAW5xd40b5ZCKE4pIFir1oG4j3DfqBdwJXE7HqrrthPtHjyW9mLu/BBxSpn0yMLmTnxmZ9PwiIlI/Kt2DWkN0f8jM5gO/B7Z293eitgOA57STuYiI1ELSvfgWAQ58O9b2a+DPZrZb1aMSEZGGlzRB3QQ8xYZ7UgC7AM8QJjyIiIhUVdIEdSBwmbuvXxDj7isIuzz8Sy0CExGRxpY0Qb0DfLRM++6A9rwREZGqS7oO6nbgVjPbCZhJmMG3L/Bd4I4axSYiIg0saYK6Mjr2u8CgqK0VuB64pgZxiYhIg0uUoNx9HfA94HtmNhBYVaGGk4iISLclLrdhZh8nlMzoHT1uAgrAfu5+Wm3CExGRRpUoQZnZxYSdJJYSNm/9J/C+6OmHahOaZEHl3UUkL5LO4jsNuMDdtwHeImzy+iFC+fcZNYpNUtY6bTpzJ0ykbdFiaG+nbdFi5k6YSOu06VmHJiINKGmC2h64N/r+eeBgd18IfBMYXYvAJH0LJk1mXVtbh7Z1bW0smFR2m0QRkZrqylZHA6LvXwX2ib7/K/DBagcl2Whb/HaX2kVEailpgnqAsA7qY4S6TieY2Qjg68BfahWcpKswcECX2kVEailpgjqfcK/po4Ty748TymycED0ndWDI6BZ6FQod2noVCgwZ3ZJRRCLSyJJOMx8FXOzuxbGeE83sLGClym3UD1XPFZE8SZqgbgT+CKy/GRHfOFbqh6rnikheJB3iewb4Qi0DERERiUvag1oHXGVm3wHmAyviT7r7gdUOTEREGlvSBPVM9CV1SLtHiEgedZqgzKxXtEks7n5pZ8dJz1bcPaK4QLe4ewSgJCUimap0D2q1mQ2ON5jZcDMrdPYD0vNo9wgRyatKCaqpTNuDhD34pE5o9wgRyauks/iKyiUt6cG0e4SI5FVXE5TUGe0eISJ5lbhgodQn7R4hInm1qQR1opnFd4zYAjjezBbHD3L3m6semaRGu0eISB5VSlALgDNK2hYCXytpaweUoEREpKo6TVDuvnOKcYiIiHSgSRIiIpJLSlAiIpJLSlAiIpJLqU4zN7N9gInA3sA84CR3n1HmuP2A66Pj3gNuBy539/YUwxURkQyl1oMys77AA8AvgG2BK4GpZrZNyXH9gd8AvwQGAJ8GTgROTStWERHJXppDfIcCfdx9nLuvdvcpwMvAcSXH7QQ87e7j3X2tu78G3A8ckmKsIiKSsTSH+IYBs0vaXgH2ije4uxOr3hv1vI4Ebq11gCIikh9p9qCageUlbcuB/p39QFTa4+fRcRNrF5qIiORNmj2oZUC/krb+wNIyx2Jm2wP3EsrN/5u7ryh3nIiI1Kc0e1CzACtpGxq1d2Bmw4AZwBxCcnq39uGJiEiepNmDegJoMrOxwHhgFGEa+X3xg8xsO2AqMMXdz08xPhERyZHUelDuvoow2WEU8A5wMTDS3ReZWUts1/TRhKq9Z5jZ0tjXz9OKVUREspfqQl13f4ky08XdfTIwOfr+RuDGNOMSEZH80VZHIiKSS0pQIiKSS0pQIiKSS0pQIiKSS6lOkpDstE6bzoJJk2lb/DaFgQMYMrqFwSOGZx2WiEinlKAaQOu06bw27iZYtw6AtkWLw2NQkhKR3NIQXwOYe/Mt65PTeuvWhXYRkZxSgmoA61au7FK7iEgeKEGJiEguKUGJiEguKUGJiEguKUE1gMKggV1qFxHJAyWoBjBkdAu9CoUObb0KBYaMbskoIhGRTdM6qAZQXOukhboi0pMoQTWIwSOGKyGJSI+iIT4REcklJSgREcklJSgREcklJSgREcklJSgREcklzeKrE6r3JCL1RgmqDrROm87cCRNZ19YGhHpPcydMBFTvSUR6Lg3x1YEFkyavT05F69raWDBpckYRiYh0n3pQPdj6Yb1Fi8s+37b47ZQjEhGpHiWoHqp0WK+cwsABKUYkIlJdGuLrocoN68VpM1gR6enUg+oh5ky8lb8/8iisWwe9eoX/dqIwaKBm8YlIj6cE1QPMmXgrf3/4kQ0Nm0hO+99+SwpRiYjUlhJUDpWuaUo62UHDeiJST5SgcqR12nTm33Yna5YsWd/W2Qy9osKggVqcKyJ1SQkqRc+edS4r33hz/eMtd9qR/cbfACSblbeRXr00nCcidUsJqgo6rEeKJjCUTlQoTU4AK994k2fPOpf9xt+wyVl55XzgiMOr9juIiORNqgnKzPYBJgJ7A/OAk9x9RpnjhgB3AJ8AWoGz3f2hNGKslGzem/1Kh5l0HzjicLbZY2jHnk80gaF0u6HS5FRUbN/UfaamQoH21as7XPt/nf7vVfqtRUTyJ7UEZWZ9gQeAccBwYBQw1cw+7O7vlRw+BXga+DxwCHC/mX3M3efVMsaNhtliyea162+E9vYNB69bx98ffoTWx5+kvZOeT3G7oST3hQoDB3R6v6lXocBuY07X/SURaShpLtQ9FOjj7uPcfbW7TwFeBo6LH2RmuwP7A5e4+yp3fxz4FXByrQOsOMwWT07x5k0MyyWdgTdkdAu9CoWN2ntv3azkJCINKc0hvmHA7JK2V4C9yhy3wN2XlRx3YA1jA2qzd11xu6Etd9qx7DDfljvtCGzYdVwlM0REgjQTVDOwvKRtOdB/M4+rukrDbJ1qaqJX375le17xdUn7jb+h4iw+CElKCUlEJEgzQS0D+pW09QeWbuZxVTdkdEvnU72bmsoO833gs59hmz2GbnIWH9AhGYmISGVpJqhZwNiStqHAT8scN8TM+rn7ithxs2ocX8dhtoSz+Ioz6dTzERGprjQT1BNAk5mNBcYTZvHtDdwXP8jd3cxeAK40s28DnwSOAQ5OI8hKw2yDRwzX1G4RkZSkNovP3VcBRxIS0zvAxcBId19kZi1mFh/CGwXsQVgDdTtwsru/lFasIiKSvVQX6kZJ5pAy7ZOBybHHbxCSmYiINCgVLBQRkVxSghIRkVyqx81iewMsXLgw6zhERGQTYp/VvUufq8cEtQNAS4sK94mI9CA7AHPjDfWYoGYA/wq8BazNOBYREamsNyE5bVTZoqm9k01QRUREsqRJEiIikktKUCIikktKUCIikktKUCIikktKUCIikktKUCIikktKUCIikktKUCIikkv1uJNERWa2DzCRUCxxHnCSu2+0grkemdmBwIPuPjjrWNJgZocDPwA+Qqgtdo2735JtVLVnZkcBVwG7EH7vqxvh9wYws22BF4FL3P2ujMOpOTM7CbgFaIs1j3H3n2QUUlU1VA/KzPoCDwC/ALYFrgSmmtk2mQZWY2bWZGanAFOBvlnHkwYz2wm4F7iC8F5/GfgPMzsi08BqzMx2AO4BvuXuWwP/BxhnZvtmG1lqJgIfyjqIFO0LXOvuzbGvukhO0GAJCjgU6OPu49x9tbtPAV4Gjss2rJq7FDiD8GHdKHYG7nb3+9x9XdRLfhL4l0yjqjF3fwsY5O4Pm1kvYACwBliSbWS1Z2ZfBbYB/px1LCnaD3g+6yBqpdES1DBgdknbK8BeGcSSponuvh8wM+tA0uLuT7n76cXHZvZ+wibCz2UXVTrcfYmZ9ScM+0wFJrj7axmHVVNmtgvwPeCkrGNJi5n1JtyqGG1mfzOzOWZ2oZk1ZR1btTRagmoGlpe0LQf6ZxBLatz9b1nHkCUzex/wK+AZwhBvI1gJbAUcAJxkZidnHE/NRB/UPwPOd/dGKgQ3iPBH508I9xuPJYyUnJFlUNXUaJMklgH9Str6A0sziEVSYGa7E5LSLKDF3ddlHFIqot9zFTDTzG4FjgHuyDaqmvku4O7+/7MOJE1RMh4Ra3rezG4CRgE3ZxNVdTVaD2oWYCVtQ6N2qTNmNpzQa7ofONbdV2YcUs2Z2Qgze7akuQD8I4t4UvIl4Fgz+4eZ/YMwZH+zmdXFh3RnzGxPM7u0pLkvofdcFxqtB/UE0GRmY4HxhL809gbuyzQqqToz2w14ELjY3W/KOp4UPQ98yMy+DtwAHAScDHwh06hqyN2Hxh+b2fPAuAaYZv4P4Btm9iahd/xx4BzgrEyjqqKG6kG5+yrgSEJiege4GBjp7osyDUxqYQywNWFq+dLY139mHVgtufs/gc8BXyT8G78VOMXdp2UamFSdu/8VOBo4DXiPsKzicne/J9PAqkgVdUVEJJcaqgclIiI9hxKUiIjkkhKUiIjkkhKUiIjkkhKUiIjkkhKUiIjkkhKUJGZmW5nZ5Wb2qpmtMLP5ZnZtVIOnWtdorta+cWbWx8zOrMa5qsHMTjSzxSlda2xxlwEze9LMftiNc23We2JmT5vZIZt73bwws/3N7LGs42hESlCSSFQz6w/AvxFWqw8DTgeOINTUKt3jcHN9g+ptdvkV4LIqnavHiGphnQtcU6VTdvk9iTboHUr4N9OjuftM4N2onIekqNG2OpLN9wPCHzSfcvcVUdt8M5sFzAG+RnU2qKxmqYC6KTvQRd8E7nH3am2CvDmv46eBp9x9TZViyNqNwI/NbFKjbDicB0pQsklmVgBaCFVaV8Sfc/c3zOwwwKNjm4CzCfuBDSHU27rI3R+Knr+LUOKkH6E8wNvA7e5+hZmdSKjpg5m1E0oILCL0BEYSygssBG5x9yui43oDFwGnAu8HZkTXHgT8OHauw9z9yZLf6/uEvRjnEeoItQH/Dxjr7mujWJvd/djYz7wO/NDdx0c/vxfwGqGHsQq4HPgT8CNgV2A6YRf1d2LnuIjQKyE67nvuvjZ67iDgOkIhujeA26LrrYten3MImxsfTSjl3qEIpZk1AycSerZxg8zsYeCw6Pcd4+5PxH7ufML7NoCwn9/57v6HzXlPIocTalEV/01cStiSZ0vgTuBjwE/c/S4z26rS+aL34T1CuZzjgHeBrxMKMV4DDCbsu/hVd1+1Oe9LtD/nGYRCl0uB3wBnxJL8b6Pr/28ap2RL5jTEJ0nsSqhUOqPck+7+e3d/O3p4EeHD6BLCh//9wK/MbJ/Yj5wCvAnsD9wOXG5m+wG/AK4FXgB2IHxAXwccTPjwMsJfssXjia5zDnAeYbPMNwgfLk9Hbe9E5/p9J7/bUYSS8AdH5xoTXSupowgfXPsSyo1fR9iI+EzgM4REc17s+AGE5PEp4ATCh/aFAGY2GHgE+C/CB+w5UTzfjP38xwnVcfcFflomnhHAWsLvH3cC8BSwD2Fz5EeiIn+Y2WmEIcEzo/M/BDwePb857wnR7/5o9P23CB/+pxAqGhsdy0QkOd9pwNzodXmEsDnqBYSS9l8ilBM5PnZ84vfFzL4MfJ+Q9D5CSPDHRNcE1pcveZSwl6ekRD0oSWK76L//rHRQ9JfyecCV7j4lav5+1Cv4JqEXBjDX3b8bfX9F9Nfrfu7+rJktBdYUC8+Z2e8If03/KTr+GjO7BBhmZn8iKmVfrAVkZmMIHzZbR/G2b6KI3QrgzGgjYTezMwgfXvdu4jUpWgmc5+5rolo83wHGu/tTUTwPAXvGjl8LHO/ubwAvmNkV0WtzJSEZzXD3y6NjXzOzbxM+sH8QO8elFX6n/YHZ7l66yeaj7n5V9P3FZnYkIWFcTPij4px3AHUAAATmSURBVEJ3/030/FVmdiihl3V+V94T4Fkz2xXYwt09ev5swiamv45+/nhCoiuqeL6obY67Xxn9/C2EHu93iz9jZn+k4+vclfflb8CJ7v5g9PgvZjYtun7cLMImvJISJShJojjzbLuKR4WhloFs/Nf7bwl/6RaVlh9fAvTp5Jw/A44ys9HA7oShoWagd3StQcR6du6+hGj4zKy09FdZf4mSU9F7FWIp5/XYfZZiteZ5sedXRnEW/TVKTkUzgA9GMyH3BA6NEkJRL6CfmQ0oXmMTCfcDbHi/4krfk5nAntGQ4BDgtuiDv6hAGPIsp9J7ArHek5kNBD5Ix/foHTN7tQvng3Cfs6iz17kQe5z4fXH3adFMvSsIEzv2JPTkSnuobxP+jUtKNMQnScwh/M95QLknzew6M7uA0Bspp4mO/9ZWdXJMOXcSJl+sIHxgHMyG4nvF83RnS/5KsZQ7b+kfdavLHFPpJvraksfF12V1dO57CR/Qxa+9CcNOxd5rZ69x/NrlXsty113FhiTw1ZLr7kGYpVlOpfcEYvef2PD6VPqs2dT54ueJq/Q6Jz4+us/2FOEe5sOEIcNflTm0Nxu/jlJDSlCySdH4+yTgbDPbMv6cme1MGKtvc/f3CMMlB5ec4pOEyRJJrE8KZrY1MBo4wd0vcvdfED5U3wc0RbWPWgn3GYo/08/M/h4NK3a3lkzxWsVzN9P9v6B3NLN4T/QThF7cMmA2MNTd5xS/CIniEip/GMctJPQqS+1d/CYaij0QmBW9hguBD5Vc9yw2TLRI/J5Ek1YOBR6D9fWp3iQMmxbP8T5C0t3k+RL+zt01BrjG3c909zuAF6P4Sq8/kPBaSUo0xCdJXQ58HnjCzL5HGKbbB7ia8D/0rdFxPwAui6p8PkuYdXUE4UMriaXA9tF9jDeAZcAXzWw+YajoGsIHR3E453rgu9HsulcJ91T+CTwH7Ag0m9kwYN5mlHyfAZxsZiMJyeNSuv8XdG/gbjO7kDAj7qLoC2ACcI6Z3Rh9vxNwC/BANIsvyfmfBS40s97FmYGRYyxU2f0NYcbjzmxYFnA1cImZvUX4nb9CSFDFiQxdeU8OILzWxUkzEN6j75jZPMIw2xXAVoTEt3IT50vD28Bh0b+TJmAs4f7Tn0uO24dOJgpJbagHJYlE03EPIdy7uIVww/hawiy9z8Y+/McTPvCuJvwPfgxwVPHmdAL3EO5JzSIMNX2FsDh4FmH4ZyphllnxL/JrCDO6bickpR2Az0f3lR6L2p4jJNeumgTcFV33KcL0699txnniZhOSyHRCUr+e8Hri7m8Skvn+hFlzkwiz6M4re6byniB8yO5f0n4TYYr0i4Q1Sp9z99bouRuAHxLes1nA/wWOdffi79qV9+RwNszeKxpHeB1/QrgXNgf4C7DK3Vdv4nxpOJeQLGcC/01IjP9Bx555L8IMxAfLnUBqQxV1ReqMmU0A1rr7OVnHAmBmnwP+FJsFuAVhIsfR7j490+ASMrMjCH98mRbqpkc9KJH6czUwKrrXkwcnAz8zs73M7COEHts/gGeyDatLzgSuUnJKlxKUSJ1x978Qhg4vyDqWyFmE+4LTCbs57AZ8xt07m8aeK2Z2AGHSxl0Zh9JwNMQnIiK5pB6UiIjkkhKUiIjkkhKUiIjkkhKUiIjkkhKUiIjk0v8AINMhYRYLcLoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_sweep_frame(frame)\n",
    "\n",
    "decorate(xlabel='Contact number (beta/gamma)',\n",
    "         ylabel='Fraction infected')\n",
    "\n",
    "savefig('figs/chap14-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that the ratio `beta/gamma`, called the \"contact number\" is sufficient to predict the total number of infections; we don't have to know `beta` and `gamma` separately.\n",
    "\n",
    "We can see that in the previous plot: when we plot the fraction infected versus the contact number, the results fall close to a curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the book we figured out the relationship between $c$ and $s_{\\infty}$ analytically.  Now we can compute it for a range of values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_inf_array = linspace(0.0001, 0.9999, 101);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_array = log(s_inf_array) / (s_inf_array - 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`total_infected` is the change in $s$ from the beginning to the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac_infected = 1 - s_inf_array\n",
    "frac_infected_series = Series(frac_infected, index=c_array);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the analytic results and compare them to the simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap14-fig02.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxcVf3/8dfMJJk2Tbom6V66AKctO4WytwVFFBGKqCy1KNsXEVD5icqiaIvI94sKCK2WTcB+q/QrO7Kr2CJr2SktH0oXQ9ckXbO3yczvj3uTTqeTdNpmlmTez8cjj2TOnHvvJyncz5xzzxKIRqOIiIhkm2CmAxAREUlECUpERLKSEpSIiGQlJSgREclKeZkOoKM558LAkcAaoDnD4YiISPtCwEBggZk1xr7R5RIUXnJ6OdNBiIjIbjkB+HdsQVdMUGsA5syZw4ABAzIdi4iItGPt2rVMmTIF/Ht3rK6YoJoBBgwYwJAhQzIdi4iIJGenRzIaJCEiIllJCUpERLKSEpSIiGQlJSgREclKGRkk4ZwbD/zNzMraeH8YcB9wNFABXGlmz6QxRBERybC0tqCccwHn3MXAC0BBO1UfAj4A+gGXAA8550amIUQREckS6W5BTQO+DPwS+GmiCs65/YEjgJPNbCvwT+fck8BFwPXpClREpCuIRqNEohCJRIlEo953/+fm5p3LIpEozbF1dlF/WP+elPbpnpLY052gZpnZDc65Se3UGQuUm1ltTNnHwPiURibSRVTMm0/57Dk0Vq0nXNKPYVOnUDZxQsbPtafib6zNiW60kbbfT1R/h5twaxl+eSSmDjE/R1rrtJ4vyZv8DnE0R2lOov7uxu6VRfy/x45xplLvojB/+sUpBAKBDj93WhOUma1OoloRUBdXVgcUdnxEItkjPhkM/eZ5FB11LJtrGtlSt5Xm5pabpH8zTXCz2rhoMWtf/AfNTb2IFvUm0hjg9Qf/Tr9lDRSOGrVbN76aVavZtHAREfYjUro/UQJE//IePd6qIb+kpO2bfDRKc3Mk5vzEJBHvJr/DDTo+nrhYpGOEggGCLV8B73so7nUwGCAU+3PrexAKBuOOh2AgwMH7laYkOUF2riRRC8S3FwuBmgzEIrJLu9PK2NbUzOaarWyqbmRTTSOb/a/Vi5ey+oOPqc07hNrB3agLdaPuyfU0P7UHY4P6Hb1z2cJ6WLhw989VNGrnstXbYPVOq9Kk1I43VQgGgwQDLTdY70YZDAX9m2t8/Z1vwtuPTXCDTlQnyfoJr9t6vB9XMBiTBLzYQ6Hg9pt+MrGH4urson4wmJoEkmrZmKAWAcOcc93NrN4vG+2Xi2SVinnzsZl3szGST023Murqu/Hag38n/HE92/oNYHNNI5uqtyei2oamtk/WY8RORQXRJvqW9KS4sIC8ULCNG9/21xv+/QpBogSiUQJECUajBIkQIMrg006NuYkSc1OMuVH6N9BgEFbccx/B1vNEWs8XCMCYH121WzftXd1Ed5VEJDdlXYIyM3POvQ/c5Jy7FjgWOAM4JrORSS6LRqOs39zAqooaVlXVsKqihpWVNSz/aCWbhpxFNBA3IPaTBmDFTucJBgP0LiqgV1GYXkVhevvfNz/2Vwqb6ilsbqBHcwOF/lc+EY679eGk43xr/r00VlbtVB4uLeGIyQft1u/81oPr2z7XIYN361wieyIrEpRzbgpwl5kV+UVnAXfjzYGqAi4ysz3onxBpW6KuueKjj2VlRbWXiCprWVVZ05qUGrcm2F4s1INANEKvbdUUN9XtkGAOuniqn4C8hNS7OEyPbvkJWwRvPXUnjRsTJ4PdMWzqFJbOnEWkcfu2OsFwmGFTp+zWeTr6XCJ7IiMJysz+BfSOeT0HmBPz+jPgS+mPTHJFxbz5LJl5F5XRbqwuHsVqSln98HKqntxIW4/le/YoYHBpkfdV5n3fdOdvKFxXTh6RHeqGS0s44ridu+za0lHJoOXZV0eMvOvIc4nsiaxoQYmkw8YtDVj5Rj4p38hbzy9h1eAz2RrM36FOMBph2MBerQkoNhkVF+48t7xiyhlZmVg6Kol05LlEdpcSlHRZ25qaWbRsA299vI63P67gs3XV29/M87rOem6rYVBDFYMaKxnUUEX/rRuZdOvcpK+RrYlFpCtQgpJOp71h3es21PH2x+t4e3EFH3xaSUPMc6NuBSH2H9aH/Yf1IfTUXyitWEZRc8MO597dZz6gxCKSKkpQkjWSmU9UMW/+Dl1qjZVVvDPrT6xd0sDbm/JZsWbLDvWHD+zJuNFljBvdn9HD+5Kf5422qyg62TtPzLgHDQAQyS5KUJIVEiWepTNnAeyQpMpnzyHS2Eh9MMzHRfvwUfEIVnbvD4vrgXoKu+Vx6P6lHO76M250GSW9E68RpgEAItlPCUqyQkviiRVpbKR89pzWpBGNRllSE+LtARNY0mMYEX/uUV6kif1qV3LW977B4aP7t7aSdkVdcyLZTQlKskJj1fo2yxsam/jXOyt5+pXlrBh8CgCBaIQRtas4oGY5+9WU07OkN0ccODCdIYtIiilBSVYIl/TbadWC+mAB7wwax503vkBt/TYAeoYDHFzxIYdsWExxs7cSlp4diXRNSlCSFWInqtYHw7zZewxv9x7jzVOq34bbpw+nHTeC4w4ZxMZXelM+u5zGqgY9OxLpwpSgJC12NUKvbOIEGpuizH10Aa8VDG+dQHvo/qWc+wXH2BH9dqirhCTS9SlBScrtaoReJBJl3rsrefCVBtZ32w+Aw10Z537BMXp434zFLSKZpQQlKdfeCL3qUQfz+0fe55PyTQDsO6QXF59xEAeM7JfoVCKSQ5SgJOUSjdDbFgjxr8gw3vzdfCKRKH17hjn/1LGcOG6o9v8REUAJStIgr6iIpurt6+CtCpfwdP/j2FDQi0A0yldOGMk3vziawm757ZxFRHKNEpSkXNTfwCJCgNf7HMjLfQ8hGghS0rSFn1z1ZT1nEpGElKAk5ZpraqkJdefJ/idQXjgAgPEbP2LChvcYPXxqhqMTkWylBCUpV9F/FP/X7VBq8grp0VTPaev+zYj6NXu0criI5A4lKNlr7c1xev71//Bg8bE0R2Fo/Vomr51Pj+YGrf4gIrukBCV7pa05TtFolOfrSvnrP5YA8LkRYY55712aIo2ES0u0+oOI7JISlOyVRHOctjVu445HPuKDgsEEgwG+e9YhnHL0PsAXMxOkiHRKSlCyV+LnODUFgjw2YCJLCwbTrSDET84/kiPG9M9QdCLSmSW3cY5IG8Il21d82BYI8fDAk1jaYyjdI1v51XePU3ISkT2mBCV7ZdjUKQTDYZoCQR4ZeCIrCgfRo7mea04qYb+hfTIdnoh0Yurik71SNnECTZEov33MWJHfnx6RRq47uT8Hf/nETIcmIp2cEpTslWg0yl/WFPFJfn+KC/O5+bsnss/AnpkOS0S6AHXxyV750zOL+dfbK+lWEGLafx2j5CQiHUYJSvbY86//h4f/uYRgMMA13zpSz5xEpEMpQcke+WjZemY9+j4Al3/tEMaN1mg9EelYSlCy29ZvrufmB9+kqTnK6RNG8oWj9sl0SCLSBWmQhLQrfp29QVPO47YPA2yu2cqh+5Vy4WkHZDpEEemilKCkTYnW2fvjn1/n456jKenVjau/OY5QSI1wEUmNtCYo59whwCzgYGAZcKGZLUhQzwF/AA4HqoFZZnZTOmOVndfZW1Y4iDd6jiYQjfCT84+kV1E4g9GJSFeXto+/zrkC4AlgLtAbuAl4wTmXaFzyHODvQF/gJOB7zrnT0xWreGLX2asLhnm67DgATtjwvnbBFZGUS2f/zCQg38xuN7NtZvYQ8BFwdoK6zv8eAKL+V0NaopRWeUVFgPfHf77saGrzujO0fi3HN63IaFwikhvSmaDGAovjyj4GDkpQ90ZgGtAILAH+bGYvpDY8iRclCsDHRcOxon0oiGzjy+teIeCXi4ikUjoTVBFQF1dWBxQmqBsFfugfcyjwVefcRakNT+I119RSHwzzYsmRAJxU9Ra9m2pprqnNcGQikgvSOUiiFugeV1YI1MQWOOeOAK4ys0F+0fvOuVuA7wL3pTxKaRUu6cfTgf2py+vOsLq1HLJlSWu5iEiqpTNBLQKuiisbDfwprmwoUOCcC5hZS19SE7AtxfHltPj5TsOmTqHh1G/w4Ss1hCLNfLHyNQJAMBxm2NQpmQ5XRHJAOhPUS0DAOXcVMAM4C2+4+WNx9V4BQsA059x0YARwNTAzjbHmlETznZbMvIs5Y88F4Pity+jbVEO4tIRhU6dQNnFCJsMVkRyRtgRlZludc1/Cmwc1HVgBTDazSufcFOAuMysyswq/3q+B7wEbgbuBO9MVa66Jn+8E8Ha3fVhZ3Uz/voVc+eOrCOeHMhSdiOSqtE7UNbOFwPEJyufgzX1qef06cEIaQ8tpsfOdABqCBbzc91AALjr9ACUnEckIrVMjOw16eLXPQTSEwgxvWs/RBw7MUFQikuuUoIRhU6cQDHvLFm3O68HbvUYDMPWk4QQCgUyGJiI5TIvFSuugh/LZc3glsD/NwRBHDc5n/BknZTgyEcllSlACeElq29jDWXjLPwkCF56vkXoiklnq4pNWD71oRCJRPnfEUAaVFGU6HBHJcW22oJxzSX+ENrP5HROOZMrqyhrmv7OSUDDA2Se7XR8gIpJi7XXx/SvudRRvdfEI0Azk+z9vJfF6epLlYlePeG7oJCIFQzn5yKH076t/ThHJvPa6+Ipjvi4APgSOAcJmFsZbxHUBOy9fJJ1Ay+oRjZVVbAl254P8QQSiEU4s2pTp0EREgHYSlJnVtnzhbX9xiZm9YWbN/vsfAJfjrQohnUzs6hFv9R5DJBDC1ZTT+NhDGY5MRMST7Ci+nnjr48Ur3o1zSBZpWT2iIZjPe732A+DoTQtp3Loxk2GJiLRKNrk8DNzvL/T6Lt6zqKOAW9l5NXLpBMIl/WisrOKDnvuyNVjAsLq1DGjcQLi0JNOhiYgAySeoK4E/AE/EHLMNuBf4UQrikhRpHRhRWUWEQOuqEUduWqStNEQkqySVoMysHvi2c+5KoGUM8sdmVtPOYZJl4rfVWFo4mM35xfTeVs3YHg0Mn/odbaUhIlkj6Ym6zrliYKr/tQKY4JwbkaK4JAXit9V4p5f3WWN800rG33uXkpOIZJWkEpRzbn/A8PZnugxv0MS5eNuxH5u68KQjxW6rsTGviOU9BpMXaWLM2g8yGJWISGLJtqB+B/yfmY0GGgHMbCowG29jQekEYrfVeL/X/gCMrllB777FmQpJRKRNySaoY4DfJyi/FTik48KRVGrZVqOZAB8WjwLg8PoVGhghIlkp2VF89UAZ8Elc+X7Alg6NSFKm5RnTC395kdq87pQ0VzPh4q/p2ZOIZKVkE9SDwEzn3Hf816XOuQPxuv7mtH2YZJuyiRNYsTQMH63l9MlH03/SvpkOSUQkoWQT1E/xFov9JxAGXgOa8Lr9rk9NaJIKm2saeWvxOoLBAJPGDcl0OCIibUrqGZSZNZnZtUAf4CDgMKCPmf0A6J3C+KSDzXt3Jc2RKIe7MvoUd8t0OCIibUqqBeWcawYGmFkl8FFM+TBgEaDd7TqJee+sBOCkcUMzHImISPva27DwXOBM/2UAuNc51xhXbR9gQ4pikw62urKGT8o30T0c4sgD+mc6HBGRdrXXxfciUAPU+q/r/Z9bvmqAN4DJqQxQOs7L760C4OgDB9KtQIvQi0h2a/MuZWZVwIUAzrkVwK/NrK7lfedcnpk1pTpA6Tjz/QQ14TANjhCR7JfsRN3bgFnOuetiypY55+5zzml/8E6gfO0WytdWU1yYz6H7l2Y6HBGRXUo2Qd0JHAA8F1M2FTgY+G1HByUd75X3VwNe915eKOk1gkVEMibZO9VpwAVm9k5LgZnNAy4FzkpFYNKxXv1wDQDHHjwow5GIiCRndz5KtzVppqAjApHUWV1Vw4o1Wyjslsch+6l7T0Q6h2QT1FN4Sx0d2FLgnBuD1/X3TCoCk47z+odrAThyzADy89S9JyKdQ7Jjja8CHgc+8OdCRfGWPHoBbzt4yWJvfOR17x190IAMRyIikrxkt3zfCEx0zo0FxgJbgU/M7ONUBid7rmLefMpnz2HjhhoWD/86oWCAw11ZpsMSEUla0rM1nXNBYDgwFLgf2N8519PMkt5uwzl3CDALb/TfMuBCM1uQoF4xXvfh6XittYeBK8xsW7LXymUV8+azdOYsIo2NLC0eSTQQYFj9WmreeI1Cba0hIp1Eslu+DwTeAx7C20G3L3ANsNh/FpXMOQqAJ4C5eAvM3gS84JzrmaD6H/06w4ExwBHAj5K5jkD57DlEGr1VqZYWepNyR1WXUz5bO6OISOexO1u+LwJK8ZY8Avgm8Lr/XjImAflmdruZbTOzh/AWnj07tpKfDM8ALjGzLWZW4b/W3TVJjZVVAEQIsLzQG1a+b+3K1nIRkc4g2QR1EjDdzFoXi/WXPboeOCrJc4wFFseVfYy3fUesw4ByYIpzbplz7jPgcmBVkteRQACAVd1KaQwV0HfrZno31bSWi4h0Bsk+gwqQeB5UCd6AiWQUAXVxZXVA/FJJffG69g7Ee1ZVhjfMvRr4VZLXym3RKADL/NbTyLpVO5SLiHQGybagHgF+63e/RQGccwcDM/GeKyWjFugeV1aItyp6rEYgBPzQzGrMbBlwK/DVJK8jvuWtCWp1hiMREdl9ySao/wesxetmK8J7HvUusNx/LxmLABdXNtovj9UydD12p17tDbEbQsVF1AXDrA33IxRpZmj9utZyEZHOos0E5Zz7L+dcEYDfkjkX2Bf4Ct7AhjFmNnk3hpm/BAScc1c55/Kdc+fgdeE9FlvJzD4E3gJuc871cM7tgzdR+C+7+8vlqpGXXMSKosEQCDCkoYL8aDOEQoy85KJMhyYikrT2Wia3A88CNTFbvi/Dm7+028xsq3PuS3jzoKYDK4DJZlbpnJsC3GVmLR/xT8UbHbgML4n+0Y9HklA2cQIb36uF8q2MqF9DuLSEYVOnUKY5UCLSibSXoFbj7QH1Jt4giR855+KfFwFgZtOTuZiZLQSOT1A+h5hh5GZWCZyXzDklMavJA7Zy+o0/YL+hfTIdjojIbmsvQX0buAH4Mt7AiJOBRDvoRvFaRJIl1q6vpWJDHUXd8xk5uPeuDxARyULtbfn+b+ALAM655cDJ/jbwkuU+/NT7ZzpwVD9CQc19EpHOKdnFYkcAOOfy/WMCce/Hz2+SDPpwqZegDt5Xez+JSOeVVIJyzh0F3I03eTZWAK+LL9TBcckeikajfLh0PeC1oEREOqtk5xf9DtgMTAaSXr1c0m/dhjqqNtVTXJjPPgMSrcMrItI5JJugDgKO9ucoSRZb6Leexo7oR1DPn0SkE0t2JYnFwKBUBiIdY9Fyde+JSNeQbAvqTuAe59ydwBLiFog1s2c6OjBJXsvuuY1V63ln+JkQKmLsCCUoEenckk1Q9/vf/yfBexokkUGxu+fWBcOsDxWRH2miaOkHMGxipsMTEdljyQ4zT7YrUNIsdvfcld3LABjYUMXqOS8w6EQlKBHpvNpMUM65wpb5Tc65+D2bdqB5UJnTWLW+9edV3bx5T0MaKmjcuL6tQ0REOoX2WkbVzrky/+cavA0D479ayiVDwiXbnzW1JKhBDZU7lIuIdEbtdfGdBGzwfz4xDbHIHhg2dQpLZ85iW+NW1oa9pDQkuoVhUy/McGQiInunvbX45iX6WbJLyxYar//5bzQF8+jXXMNBl12orTVEpNPTTrVdQNnECYRCQ+CxDzl4/GjKJo7LdEgiIntNCaqTa5kD9e+gg+JRDNqqwREi0jVo+Hgn1jIHqrGyijX+86f8vz9Jxbz5GY5MRGTvKUF1Yi1zoBqC+awv6E0o2kxp9TrKZ8/Z9cEiIlku2e02SoGfAuOAfHbeD2p8x4cmu9IyB2pduC8AZY0bCRHZYW6UiEhnlewzqPuAo4D/RdttZI1wST8aK6tYGy4BvBUkWspFRDq7ZBPUBOAMDTfPLi1zoFqeP/Vv3EAwHGbY1CkZjkxEZO8lm6A24m1YKFmkZa7TusdWAjCssJlRl3xHc6BEpEtINkHdAMxwzl1F4u02tBZfhhSOP4aNTz5Dfl6QL/7hFvJCGvciIl1DsgnqVqA38Hob72u7jQxZvspr2A4f2FPJSUS6lGQT1NdSGoXssaV+gho5uFeGIxER6VjJ7gc1D8A51x3YD2/+1FIz00rmGbZ8tZegRilBiUgXk+w8qBDwK+D7bJ8HtdU59wBwhZk1pSxCadcyvwU1QglKRLqYZB9a3AR8EzgfGOp/nQ+cijeAQjJgW1Mzn62rJhCA4QN6ZjocEZEOlewzqPOBi83smZiy/3POVQN3oySVEeVrq2mORBlc2oNuYa37KyJdS7ItqCLg0wTly4CSjgtHdseKNd6iHsMHqXtPRLqeZBPUAuDyBOVXAG93XDiyO1oS1IiB6t4Tka4n2X6hnwD/cs5NYvtcqKOB4cAXk72Yc+4QYBZwMF7r60IzW9BO/Xz/ek+Z2S+SvU5X17IH1Hv5h0HhQPps+AxwmQ5LRKRDJdWCMrO3gMOBF/EGSJQATwGjzey1ZM7hnCsAngDm4k36vQl4wTnX3sf/XwKHJnP+XFExbz6f3jGTxsoqKgt6A9D06J+1B5SIdDlJP1k3s0+Aq/fiWpOAfDO73X/9kHPuCuBs4J74yn5r7WTg+b24Zpez/J4/Em1qoi4YpjavOwWRbfRs2Mzye/6oNfhEpEtpM0E5594ETjGzjc65BUC0rbpJ7gc1FlgcV/YxcFCCa/fBS1pn4rWixNdU7c2Nrgx7raeSrZsIxJSLiHQV7bWgngYa/Z//1gHXKgLiF5WtAwoT1J0F/N7MFjqnZyuJVBb0AaC0cVOGIxERSY02E5SZTYt5+RLwmplti63jnAvjTdZNRi3QPa6sEKiJO+e38Z5x3Y7sJFRcRHN1DVUF3tDykq2bWstFRLqSZIeZv4Q3sCHeSODPSZ5jETsPNRvtl8c6FxgPbHTObQK+DFzjnOuIVlynN/KSiyAUoqpgexcfoZBXLiLShbT3DOoyoKUVFQAWOefin0MVAe8mea2XgIC/p9QM4Cy84eaPxVYys1Pi4ngceE/DzD1lEycQjUZZ/0QFAIOKg+z3zSs0QEJEupz2nkHdg9ctFwT+CNzIjrvqRvG65/6RzIXMbKtz7kt4z5emAyuAyWZW6ZybAtxlZuqnSkLBuKOof+p5enTL48Tf/I5AIJDpkEREOlx7z6CagD8BOOeWA68CxWa2wS87Enh3d1YyN7OFwPEJyucAc9o4ZnKy588Vn63zRuwN6V+s5CQiXVayz6AqAQOujSl7CvjQOTeqw6OSdn22zhtXMqx/cYYjERFJnWQT1J3Ay2x/JgUwAngD73mSpNHKlhZUmRKUiHRdySao8cB0M2sdEm5m9XiTaI9LRWDStpUV3j/D0P56ZCciXVeyCWoDcGCC8v0BLWGQZisrvD/54DIlKBHpupJdi+9e4G7n3FDgLbwRfIcDPwPuS1FskkBdwzaqNjeQFwrSv2+PTIcjIpIyySaom/y6PwNK/bIK4Dbg1ymIS9qwurIWgIElPQgFNYJPRLqupBKUmUWAnwM/d86VAFvNbEtKI5OEVlV6z5+GqHtPRLq4pLfbcM4dhrciech/HQDCwDgzuzQ14Um81X6CGlyqBCUiXVtSCco5dz3eShI1QA+8FSV6+W8/k5rQJJGVrQlKz59EpGtLdhTfpcCPzKwnsAZvDb3BeNuxt7llu3ScinnzeeviS/nktQ8AKFy1NMMRiYikVrIJagDwiP/ze8AxZrYW+DEwNRWByXYV8+azdOYsGiqr2JjvTc6tmztb27yLSJe2O0sd9fN//gQ4xP95FTCoo4OSHZXPnkOksZG6UDcaQwWEmxvpVr+F8tkJly8UEekSkk1QT+DNgzoUb9uM851zE4H/B/wnVcGJp7FqPUBr66nPtmoCMeUiIl1RsgnqarxnTQfibf/+T7xtNs7335MUCpd4jdcN+T0B6LuteodyEZGuKNkEdRZwvZn9r5lFzezbeDvslpiZRvGl2LCpUwiGw2zyW1C9t20hGA4zbOqUDEcmIpI6yc6DugN4E2jtU4pdOFZSq2W33Cce/sh73S3KqG9/R7voikiXlmwL6g3gzFQGIu0rmziBxlFjATj2x1cqOYlIl5dsCyoC/Mo591NgOVAf+6aZje/owGRna6u8dfgGlBRmOBIRkdRLNkG94X9JmlXMm0/57Dls3lBNzYizCYegd1E402GJiKRcmwnKORf0F4nFzKa1VU9Sp2WCbqSxkY1hb8Rez/pNVM5/WV18ItLltfcMaptzriy2wDk3wTmnj+9p0jJBF2BTvrc4bO+tmqArIrmhvQSVaLOhv+GtwSdpEDsRd1NeyxDzak3QFZGckOwovhbaIS+NYifitragmmo0QVdEcsLuJihJo5YJugCb/QTVh0ZN0BWRnJD0hoWSfi0DIcpnz2ldReKgb5ymARIikhN2laC+7ZyLXTEiD/imc64qtpKZ/b7DIxPAS1L9TjiB6p88BZEoY05RchKR3NBegioHLosrWwtcEFcWBZSgUmjD5gaaI1F6F4cJ54cyHY6ISFq0maDMbHga45B2rNvgrSDRv69WkBCR3KFBEp1AxcY6APr3UYISkdyhBNUJrNvgLX3Yv58SlIjkDiWoTqDSb0GVqgUlIjkkrcPMnXOHALOAg4FlwIVmtiBBvXHAbX69LcC9wI1mFk1juFmjpYuvrE/3DEciIpI+aWtBOecKgCeAuXi78d4EvOCc6xlXrxB4Gvg/oB/wOeDbwCXpijXbVGz0uvjK1IISkRySzi6+SUC+md1uZtvM7CHgI+DsuHpDgdfMbIaZNZvZEuBx4Pg0xpo1IpEolX6CKu2tFpSI5I50dvGNBRbHlX0MHBRbYGZGzO69fsvrS8DdqQ4wG22qaaSpOUJxYQHdwlr4Q0RyRzpbUEVAXVxZHdBmv5W/tcdf/HqzUhda9moZIFHWV60nEckt6fxIXgvE32ULgZoEdXHODQAewdtu/vNmVp+oXldXucn7tcsGnyYAABQ6SURBVEt6KUGJSG5JZwtqEeDiykb75Ttwzo0FFgCf4iWnjakPLzu1PH8q0yoSIpJj0tmCegkIOOeuAmYAZ+ENI38stpJzrg/wAvCQmV2dxviyklpQIpKr0taCMrOteIMdzgI2ANcDk82s0jk3JWbV9Kl4u/Ze5pyrifn6S7pizSZVmzSCT0RyU1qHhZnZQhIMFzezOcAc/+c7gDvSGVc2a2lBlWqSrojkGC11lOVaWlAlakGJSI5Rgspi25qa2VTdSDAYoE/PbpkOR0QkrZSgstj6zQ0A9C0OEwoGMhyNiEh6KUFlsZYEpe49EclFWjsny1TMm0/57Dk0Vq3nk4EHQuFh9NMQcxHJQUpQWaRi3nyW3H4nRCIAbKyLQCEU1qzPcGQiIumnLr4ssvT3d7UmJ4DqPK/lFHnnjUyFJCKSMUpQWSTS0LDD6+q8HgAUNWzJRDgiIhmlBJXFqvO89feKmuIXgRcR6fqUoLJYdchLUMVKUCKSg5SgslSEADX+M6iiZiUoEck9SlBZJFxa0vpzXShMNBCke3MDPUr6ZjAqEZHMUILKIsOmTiEYDgNQ0/L8qbmBYVOnZDIsEZGM0DyoLFI2cQIA5bPnUF3nde/1H1LSWi4ikkuUoLJM2cQJlE2cQOVrK+Dh9xk4fFCmQxIRyQh18WWpDS0LxfbSKuYikpuUoLLUxmo/QWmbDRHJUUpQWap1qw0lKBHJUUpQWWrDFiUoEcltSlBZapPfxdenWAlKRHKTRvFlUOzeT+GSfgybOoWyiRNojkTZVN0IQO/icIajFBHJDCWoDKmYN5+lM2cRafQSUWNlFUtnzgIg/7DxRKLQs0cB+Xlq5IpIbtLdL0OW3/PH1uTUItLYSPnsOWz0W096/iSS22644QZGjx7NkiVLOuR811xzDTfddNMeH//kk09yzjnndEgsyVCCyoCKefNpqq5O+F5j1frWARJ91L0nkrNqa2t59tlnOfPMM5k9e3amwwHg9NNP56GHHkrb9ZSgMqB89pw23wuX9GOjn6D0/Ekkd/3tb39jzJgxXHzxxTz11FNs3rwZgEcffZQLLriAa6+9lnHjxvH5z39+h6Tx5ptvct5553HMMcdw2GGHcfnll1Md94F47dq1jBkzhvLy8tayxx9/nK9+9asA3H333UyYMIGjjjqKKVOm8MEHH7Re+4wzzgBgw4YNXHrppRx55JFMmjSJa6+9loa4TVf3lp5BZUBj1fo23xs2dQrL1MUnklLT7n2dtxavS9v1jhjTn59ffPRuHTN37lwuuOACRo0axYEHHsjDDz/MRRddBMCrr77KtGnTuPHGG3nssceYPn06p556Knl5eVx++eVMmzaNU089lYqKCr71rW8xd+5cLr744tZzDxgwgPHjx/P0009z2WWXAV5CPP3001m4cCH33XcfTzzxBGVlZdx5553ceuutPPDAAzvEN3PmTIqLi3n11Vepra3l/PPP57nnnmPy5Ml798eKoRZUGn06625eOfPrEI0mfD9UXETZxAlsqmkZwacEJZKLPvzwQ9asWcMpp5wCwLnnnsucOXOIRCIAlJaWcs4555CXl8fkyZPZunUra9asIRwO8/DDD3PqqadSV1dHZWUlffv2paKiYqdrnH766TzzzDOA1xp68803Oe200+jRowe1tbU8+uijLF26lCuvvHKn5ARQVFTEwoULeeGFF4hGozz++OMdmpxALai0+XTW3ax79vk23w+Gw4y8xPt0tFHPoERSandbM+k2d+5cqqurOfHEEwGIRCJs2LCBf/7znwD069evtW5+fn5rnVAoxPz587n//vuJRCKMHj2aLVu2EE3wofiUU05h+vTpLFmyhDfffJPx48dTUlJCSUkJM2fO5P777+cPf/gD/fr14/LLL+frX//6DsdffvnlBINBZsyYwdVXX824ceOYPn06I0eO7LC/gxJUirTOcaqsgmAQ/E8+iYRLS1rnQAGto/j0DEok99TU1PD0009zzz33sO+++7aWz5o1i9mzZ7c+A0rk3Xff5fbbb+evf/1ra6Jo6cKLV1RUxEknncTzzz/PggULOOusswCoqKigb9++PPDAA9TX1/Pcc89xzTXXcPzxx+9w/CeffMI555zD97//fdasWcPNN9/M9OnTE7a29pS6+FKgZY5TY2WVV9BOcgI44t67dtjzqWWSrlpQIrnnySefZMCAARxzzDGUlpa2fp199tm8/vrr1NfXt3lsdXU1wWCQcDhMJBLh2Wef5eWXX2bbtm0J659xxhk8++yzLF68mM9//vMALF26lIsvvphPPvmE7t27069fPwoKCujevfsOxz744IPcdNNN1NbW0q9fP7p160avXr067g+BWlAd4u0rvk/DZyu3F4RC0Nyc3MHBnT8jtCxzpGdQIrln7ty5nHbaaTuV77///hxwwAHcfPPNjBo1KuGxJ5xwAl/5yleYPHkywWCQMWPG8I1vfAMzS1j/+OOP57rrrmPSpEkUFnq7eB9zzDFceumlXHrppWzcuJFBgwZx22230bt37x2Ovfbaa7nhhhs48cQTaWpqYvz48UybNm0vf/sdBRL1TaaKc+4QYBZwMLAMuNDMFiSoNwy4DzgaqACuNLNnkrzGcGD5P/7xD4YMGbLXMX86627WPf+i1woKBAgUFBBtbNzebbeL7rtd6f+lU9j3O//V+nrrtmbOuuZvhIIBHrvlKwQCgb3+HURE2nL66adzzTXXcOyxx2bk+itXruRzn/scwAgzWxH7XtpaUM65AuAJ4HZgAnAW8IJzbh8z2xJX/SHgNeDLwPHA4865Q81sWarjrJg3n+X3/DHxRNpo1EtOsD0p7WlyCgbpf8rJOyQnoHUEX6+isJKTiKRMeXk5r7/+OrW1tRx9dHYOGklnF98kIN/MbvdfP+ScuwI4G7inpZJzbn/gCOBkM9sK/NM59yRwEXB9KgOsmDefT++YSbSpKWXXCIbDjLr8Ozs8c4qlRWJFJB1uueUW3n33XX7zm98QTPCoIRukM0GNBRbHlX0MHJSgXrmZ1cbVG5/C2ABvhYeOTE7h0pIdRvHFj9ZLZHONEpSIpN6MGTMyHcIupTNBFQF1cWV1QOEe1utw7a3wsLu6DR3CuBm/2+3jhvYvpqxvIcceNLDDYhER6YzSmaBqge5xZYVAzR7W63Dhkn7bh4bvhT1NTgAD+vXgvutP3usYREQ6u3QmqEXAVXFlo4E/Jag3zDnX3czqY+otSnF8DJs6pf1nUAlG8SXTbSciIrsvnQnqJSDgnLsKmIE3iu9g4LHYSmZmzrn3gZucc9cCxwJnAMekOsCWJBM7ii9UXMTISy5SAhIRSbO0JSgz2+qc+xLePKjpwApgsplVOuemAHeZWZFf/Szgbrw5UFXARWa2MB1xlk2coGQkIpIF0rqShJ9kjk9QPgeYE/P6M+BLaQxNRESyTHYOfhcRkZynBCUiIllJCUpERLJSV1zNPASwdu3aTMchIiK7EHOvDsW/1xUT1ECAKVOmZDoOERFJ3kBgaWxBV0xQC4ATgDVAkpsyiYhIhoTwktNOWy+ldT8oERGRZGmQhIiIZCUlKBERyUpKUCIikpWUoEREJCspQYmISFZSghIRkaykBCUiIllJCUpERLJSV1xJImnOuUPwNlA8GFgGXGhmO81mzmXOuZOB/wb2w9tA8tdmdldmo8o+zrnewAfADWb2QIbDyTrOuYHAH4ATgQbgbjP7WWajyi7OuaOBOwAHVAL/bWb3ZjaqzMrZFpRzrgB4ApgL9AZuAl5wzvXMaGBZxDk3FHgE+CXe3+hc4Gbn3CkZDSw7zQIGZzqILPYE3vJj/YGjgW85587LbEjZwzkXxPsb3WFmvfD+X5vhf4jOWbncgpoE5JvZ7f7rh5xzVwBnA/dkLKrsMhz4s5k95r9e4Jz7F3Ac8Hymgso2zrlvAT2BDzMdSzZyzh0FjASOM7NtwHLn3CSgPqOBZZc+QBkQcM4FgCjQBGzNaFQZlrMtKGAssDiu7GPgoAzEkpXM7GUz+07La+dcX7yFeN/NXFTZxTk3Avg5cGGmY8li4/CS9y+cc6ucc0uBM81sTYbjyhpmth6YATwIbMNbOPU6M4u/R+WUXE5QRUBdXFkdUJiBWLKec64X8CTwBl5XRM5zzoWA/wWuNjNtQNa2lg822/BaUl8FrlYX33Z+F18DcB7QHa+H5+fOuS9kMq5My+Uuvlq8/xBiFQI1GYglqznn9sdLSouAKWYWyXBI2eJngJnZo5kOJMs1AlvM7Bf+6/edc/fiJao/Zyyq7PJVvC7QH/mv5znn7gMuBV7IXFiZlcstqEV4o2VijfbLxeecm4DXanoc+JqZNWQ4pGxyDvA159wm59wmvO7h3zvnfp/huLLNx0ChPzCpRS5/OE5kKBCOK2vCa3XmrFz+j+QlvAeSV+H1/Z6FN9z8sXaPyiHOuVHA34DrzezOTMeTbcxsdOxr59x7wO0aZr6TF/GGTf/WOfdDvA+GFwGXZTSq7PIC3gjZ/8IbpHU4cAlwcUajyrCcbUGZ2VbgS3iJaQNwPTDZzCozGlh2uRwoxvsfpybm638yHZh0Hn6reyLe86c1wHPALWb2SEYDyyJm9hFeN9+lwCa8rs9rzCynn/dqR10REclKOduCEhGR7KYEJSIiWUkJSkREspISlIiIZCUlKBERyUpKUCIikpVyeaKudAHOuR7ANXir0A8F1gKPAjea2aYOukYRcLaZ3dcB58oHLjGzrFhtwjn3beA3ZlaShmtdBfQ2s5/vot7VQA8zm5bqmCS7qQUlnZa/d9frwOeB7+GtUP8d4BS8vb3i11rcUz+k41Y9OA+Y3kHn6jT8vcW+D/w6iep3Auf7a0BKDlOCks7sv/H+Gz7JzJ4zs+Vm9jzeCiGHABd00HUCHXSejj5XZ/Jj4GEz2+VizGbWiLftxHUpj0qymlaSkE7JORfG24L+J2Y2K8H7x+KtNL7e3wDuSuAKYBje4qXXmdkzft0H8LZa6Q58DVgP3Gtmv/S7wO6POfUIvHXlfg1MBkrxuhXvMrNf+ucL4d1cL8HbamKBf+1SvDUgW5xoZv+Ki/sXeGtCLsPbY6oR+CtwlZk1+7EWmdnXYo5ZgddNN8M//iBgCV6rbytwI/AO3pbrI4H5eKvSb2jp4gNuxWsp4tf7uZk1++c/yn9/HPAZ3lpxvzGziH/89/AWWT4dbwmjX8b9TkV4SxydYmav+mXdgV/h7RzbHfgX8F0zW+W/fyDwNjBEy4/lLrWgpLMaibeL7YJEb5rZq/4mcOAli2nADXg3/8eBJ+O2074YWAkcAdwL3OicGwfMBX4LvA8MxLtB3wocg5egHHBHTH3863wP+AFwmH/M08BrftkG/1yvtvG7nQb09q9xA96aiJOT+JvEHl+Et+DoLD/eGcB3gS/gJZofxNTvh9ctehJwPt56cNcAOOfK8HZPfg4v8X3Pj+fHMccfBlT71/tTgngmAs14v3+LWcAZwDfxtoDvhve3BsDMFuJ9EPjcbvze0sUoQUln1cf/vrm9Sn7r6QfATWb2kJl94u9L9CI73mSXmtnPzPNLvCQyzszq8fYIazKztX6r4hXgIjN708yWmdmv/Tpj/etdBvzSzB41syV4N/RH8Bbe3QxE/XO1tZ13PV5rwszsHrzkOK6Nuok0AD8ws6V4z3NCwAx/h+RXgGeAA2LqNwPfNLP3zexp4Jd4yQw/9gVmdqOZLTGzZ4FrgavjrjnNzD41s/IE8RwBLDazKLRufnkeXqvw7/6usd8B/h23Jcci/1jJURrFJ51Vlf+9T7u1oAwoYcdP7wD/Br4e83pJ3PvVQH4b5/xf4DTn3FRgf+BQvBZLyL9WKTEtOzOrxu8+cy5+C7KE/hOXvLa0E0siK8ysyf+5ZdfoZTHvN/hxtlhlZp/FvF4ADHLO9cZLZJOcc7HPjoJAd+dcv5Zr7GJH4f5s//cC72+Wx45/o+X4rbYY6/H+/SRHKUFJZ/Up3g3sSBJ08znnbsV77nFXG8cH2LEHIVFrpq0BDX/EGzn4J//rMrxWTux59ubhbnuxJDpv/P/HiTa5a28X5Oa41y1/l23+uR8BfprguJbWa3075265duzfMtm/UQgvmUqOUhefdEr+tvOzgSudc91i33PODcd7jtJoZluA1XjPc2IdizdYIhmtN1LnXDEwFTjfzK4zs7l4N9xeQMDMNuMN3jg85pjuzrl1/mCDvR2V1HKtlnMXsfetjCHOudiW6NF4rbhaYDEw2u+++9TMPgXG4D0bay/pxVqL16pssQwvKcb+jfZxzm1wzg2OqVfiHys5Si0o6cxuBL4MvOSc+zleN90hwC3AB8Ddfr3/BqY751bijQw7G29QwKQkr1MDDHDOjcQb8FALfNU5txwYhDeiL8D2LbtvA37mj677BG8zzM3Au8AQoMg5NxZY5m/mtzsWABc55ybjJY9p7NwC2l0h4M/OuWvwRilex/Yh3jOB7znn7vB/HorXKn3CH8WXzPnfBq5xzoXMrNnMqp1z9wK3Oueq8br/fgt80DKKz3cw3uAOyVFqQUmnZWYbgOOBt/BumovwbnSPA1+MufnPwEtatwAf4o0eO83MXk7yUg/jPZNahPe86Ty8Lr5FeF18L+ANPGgZyPBr4D680YDv4o3Y+7L/XOkfftm7eMl1d80GHvCv+zLwHt6gjb2xGC+JzMdL6rfhd42a2Uq8ZH4EXjfmbLzRdj9IeKbEXsJL4LEDHn7olz+GN5qxGu+DAwDOuQPwBpW8uCe/kHQNmgclIinnnJsJNJvZ95KsfzNQZmYXpTYyyWZqQYlIOtwCnOUPMW+XP4n3PP8YyWFKUCKScmb2H7yuwx8lUf0K4EEzs9RGJdlOXXwiIpKV1IISEZGspAQlIiJZSQlKRESykhKUiIhkJSUoERHJSv8f5zfy495EAAsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_sweep_frame(frame)\n",
    "plot(frac_infected_series, label='Analysis')\n",
    "\n",
    "decorate(xlabel='Contact number (c)',\n",
    "         ylabel='Fraction infected')\n",
    "\n",
    "savefig('figs/chap14-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The agreement is generally good, except for values of `c` less than 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  If we didn't know about contact numbers, we might have explored other possibilities, like the difference between `beta` and `gamma`, rather than their ratio.\n",
    "\n",
    "Write a version of `plot_sweep_frame`, called `plot_sweep_frame_difference`, that plots the fraction infected versus the difference `beta-gamma`.\n",
    "\n",
    "What do the results look like, and what does that imply? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def plot_sweep_frame_difference(frame):\n",
    "    for gamma in frame.columns:\n",
    "        column = frame[gamma]\n",
    "        for beta in column.index:\n",
    "            frac_infected = column[beta]\n",
    "            plot(beta - gamma, frac_infected, 'ro')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7xUdbnH8c8GcSehRwU5aooXqsfwVmpq5xRY1unYRVEssx1FiuE1Q6tzykt5wc4pTVNRvFuEYkdL06woUygrhcor+qigopkBagmiG4R9/nh+A8M4e/bae8+aPbPn+369eLH3b9aseWb2mvWs32X9fi0dHR2IiIjUmwF9HYCIiEg5SlAiIlKXlKBERKQuKUGJiEhd2qCvA6g2M2sF3g38DVjdx+GIiEhlA4GtgLnu3l78QL9LUERy+m1fByEiIt3yPuB3xQX9MUH9DWDGjBlsueWWfR2LiIhU8Pzzz9PW1gbp3F2sPyao1QBbbrkl22yzTV/HIiIi2byhS0aDJEREpC4pQYmISF3qkyY+M9sbuM3dh3fy+AjgKmBfYDFwgrvfXsMQRUSkj9W0BmVmLWY2EZgFbFhh05nAA8BQ4ChgppntWIMQRUSkTtS6ie8M4Bjg7M42MLO3A3sBp7v7Snf/DfBT4MjahCgiIvWg1k1809z9dDPbr8I2o4BF7v5KUdmjwN65RiYi0k8tnj2HRdNn0L70BVqHDWXE+DaGjxldd/ssVdME5e7PZdhsCLCipGwFMLj6EYmI1J9qnvwXz57DgqnTWNMekzS0L1nKgqnTAOpqn+XU4yi+V4CNSsoGA8v7IBYRkYoWz57DvImTuHvsocybOInFs+f0en8Lpk6jfclS6OhYe/Lv6X4XTZ+xNpEUrGlvZ9H0GT2OMY99llOPCWo+MMLMipPUTqlcRKRuVDuZQPVP/u1LX+hWeV/ts5y6S1Du7sD9wBQzazWz9wMHAdf1bWQi0uiqXdvJoyZR7ZN/67Ch3Srvq32WUxcJyszazKy4CW8c8A7iHqgrgSPd/aE+CU5E+oU8ajt51CSqffIfMb6NAa2t65UNaG1lxPi2Hu0vr32W0yc36rr7XcCmRb/PAGYU/f4McEDtIxOR/qpSbaenHfutw4ZGwitT3lMjxretNwABenfyL7y3ao64y2Of5fTHyWJFRN4gj9pOtZMJ5JdQqp088thnKSUoEWkKedR28qpJ1OLk3wiUoESkblXzfqA8ajugZJInJSgRqUvVvhm0Vv0m1VCLWRoagRKUiNSlPAY1NEJtp1azNDSCuhhmLiJSqlY3g9abWs3S0AiUoESkLtXqZtB606yJuRwlKBGpS7W6GbTeNGtiLkcJSkTq0vAxoxl53NG0bjEMWlpo3WIYI487ut/3wzRrYi5HgyREpG41yqCGRpyloREoQYmI9FBeI+4aITHXgpr4RER6SCPu8qUalIhURTPeXKoRd/lSDUpEei2PpSwagUbc5UsJSkR6rVmbujTiLl9q4hORXmvWpi6NuMuXEpSI9FoeS1k0Co24y4+a+ESk1xqlqWvx7DnMmziJu8ceyryJk/p9H1mjUw1KRHqtEZq6NEt441GCEpGqqPemrjyW75B8qYlPRJpCsw7kaGRKUCLSFHTPUuNRghKRptAoAzlkHfVBiUhTaISBHLI+JSgRaRr1PpBD1qcmPhERqUtKUCIiUpeUoEREpC4pQYmISF3SIAkRqVvNuAiirKMEJSJ1SXPnSU0TlJntDkwDdgMWAke4+9wy2xlwKbAHsAyY5u5TahmriPQtzZ0nNeuDMrMNgVuAG4BNgSnALDPbpMzmM4BfA5sDHwC+aGYH1ipWkf6uEZad0Nx5UstBEvsBg9z9Andf5e4zgYeBw8psa+n/FqAj/XutJlGK9HOFprP2JUuho2Nt01m9JSnNnSe1TFCjgEdKyh4Fdi2z7VnAGUA78DhwnbvPyjc8keZQqemsnmjuPKllghoCrCgpWwEMLrNtB3Byes47gUPM7Mh8wxNpDo3SdDZ8zGhGHnc0rVsMg5YWWrcYxsjjjlb/UxOp5SCJV4CNSsoGA8uLC8xsL2Cyu2+diu43s28DxwJX5R6lSD/XOmxoNO+VKa83mjuvudWyBjWfdX1LBTul8mLbAhuaWUtR2evAqhxjE2kaajqTRlHLGtSdQIuZTQYuBsYRw81/UrLd3cBA4AwzOxPYAfgyMLWGsYr0W8287IRu/G0sNUtQ7r7SzA4g7oM6E3gKGOvuS8ysDbjM3Ye4++K03XeALwIvAZcDF9UqVpH+rhmbznTjb+Op6Y267v4Q8N4y5TOIe58Kv/8ReF8NQxORfk43/jYeTRYrIk2hUUYvyjpKUCLSFHTjb+NRghKRpqDRi41Hs5mLSFNo5tGLjUoJSkSaRjOOXmxkauITEZG6pBqUiFSFboKVaus0QZlZ5iPL3etrnn4RqSndBCt5qFSDuqvk9w5ifaY1wGpgUPp5JeVnJBeRJqGbYCUPlfqgNi7693ngQeA9QKu7txLLYMwFJucdpIjUN90EK3notAbl7q8Ufjazs4Bx7j636PEHzOw44OfAZblGKSJ1rZGW8JDGkXUU3ybEDOOlNkYDLUSanm6ClTxkTS43AtekpTL+QvRF7QN8F/hBTrGJSIPI6yZYjQxsblkT1AnApcAtRc9ZBVwJfCWHuESkwVT7JliNDJRMCcrdXwUmmNkJrFsV91F3X17haSIiPaaRgZJ5Jgkz2xgYn/49BYw2sx1yiktEmpxGBkqmBGVmbwecWOH2GGLQxOHA/Wb2b/mFJyLNSstjSNYa1PeAH7n7TkA7gLuPB6YTS7OLiFSVRgZK1kES7wFOLFP+XeD+6oUjIhK0PIZkTVCvAsOBx0rK3wa8XNWIREQSLY/R3LI28X0fmGpm70m/b2FmBxJDz2fkEpmIiDS1rDWoU4nJYn8DtAJ/AF4HLgFOySc0ERFpZplqUO7+urt/DdgM2BV4F7CZu38J2DTH+EREpEllqkGZ2WpgS3dfAjxcVD4CmA8MySc8ERFpVpUWLDwcODj92gJcaWbtJZttB7yYU2wiItLEKjXx/QpYDhSW3Xg1/Vz4txy4BxibZ4AiItKcKq0HtRQ4AsDMngK+4+4rCo+b2Qbu/nreAYqISHPKOsz8fGCamX29qGyhmV1lZlruXUREqi5rgroI2Bn4RVHZeGA34LxqByUiIpI1QX0M+Ly7/7lQ4O6zgUnAuDwCExGR5pZ5uQ3gTZ2Ub1iNQERERIplTVC3ElMd7VIoMLN3EE1/t+cRmIiINLesUx1NBm4GHkj3QnUQUx7NIpaDz8TMdgemEX1XC4Ej3H1ume02JpLfgem1bgSOd/dVWV9LRGpr8ew5mnlcqirrku8vAWPMbBQwClgJPObuj2Z9ITPbELgFuAAYTfRdzTKz7dy9dEb0q4FBwPZE0+LPga8A52R9PRHpXLWTyeLZc1gwddraJdrblyxlwdRpAEpS0mPdWfJ9AJEwtgXmAJuY2SbdeK39gEHufoG7r3L3mcS0SYeVvM5WwEHAUe7+srsvTr9r1nSRKigkk/YlS6GjY20yWTx7To/3uWj6jLXJqWBNezuLputrKz2Xdcn3rYD7gJnECrqbA/8NPJL6orIYBTxSUvYoMflssXcBi4A2M1toZs8AxwF/zfg6Iv3K4tlzmDdxEnePPZR5Eyf1KpFAPsmkfekL3SoXyaI7S77PB7YgpjwC+Azwx/RYFkOAFSVlK4DSG303J2pquxB9VWOIvqivZnwdkX4jj9pOHsmkddjQbpWLZJE1QX0AONPd1152pWmPTgH2ybiPV4CNSsoGE3P6FWsHBgInu/tyd19ILC1/SMbXEek38qjt5JFMRoxvY0Br63plA1pbGTG+rcf7FMmaoFoofx/UMGLARBbzASsp2ymVFysMvCheZyrraEORfiWP2k4eyWT4mNGMPO5oWrcYBi0ttG4xjJHHHa0BEtIrWU/8NwHnmdmniWHfmNluwFRiZF4WdwItZjYZuJgYxbcb8JPijdz9QTObB5xvZuOJJDgZuCLj64j0G63DhkbzXpnyniokjWoPCR8+ZrQSklRV1gR1EpEgCgMV5hPDwG9Nj3XJ3Vea2QHEfVBnAk8BY919iZm1AZe5e2Hhw48QfVsLiVre1cTwdJGmMmJ823rDt6E6TWdKJtIIWjo6Oso+YGZfAK5z9+VFZTsC7yAS2yPu/lhNouwGM9seePKOO+5gm2226etwRHpNN8BKf/bss8+y//77A+zg7k8VP1apBnUBcYPs8qIl3xcStRoRqRHVdqRZVUpQzxFrQN1LDJL4ipmVjrgDwN3PzCM4ERFpXpUS1ATgdOCjxMCIDwHlVtDtIPqUREREqqbSku+/A/4DwMyeBD6UloEXERHJXdbJYncAMLNB6TktJY+XzhAhIiLSK5kSlJntA1xOTD9UrIVo4htY5bhERKTJZb0P6nvAP4GxQOnSGCIiIlWXNUHtCuzr7g/mGYyIiEhB1gT1CLA1oAQlImXphmKptqwJ6iLgCjO7CHickgli3f32agcmIo1DK+pKHrImqGvS//9b5jENkhBpcpWWBVGCkp7KOsw889LwItJ8tKKu5KHTBGVmgwv3N5lZ6aq369F9UCLNLY9lQUQq1YyWmdnw9PNyYFmZf4VyEWliWlFX8lCpie8DwIvp5/fXIBYRaVB5LYIoza3SXHyzy/0sIlKOlgWRass6ik9E+hHdsySNQAlKpMnoniVpFBo+LtJkKt2zJFJPlKBEmozuWZJGkXW5jS2AU4E9gUG8cT2ovasfmojkQfcsSaPI2gd1FbAP8EO03IZIQxsxvm29PijQPUtSn7ImqNHAQRpuLtL4dM+SNIqsCeolYsFCEekHdM+SNIKsCep04GIzm0z55TY0F5+IiFRV1gT1XWBT4I+dPK7lNkREpKqyJqhDc41CRESkRNb1oGYDmNlGwNuI+6cWuLtmMhcRkVxkvQ9qIHAOcCLr7oNaaWbXAse7++u5RSgiIk0p60wSU4DPAJ8Ftk3/Pgt8hBhAISIiUlVZ+6A+C0x099uLyn5kZsuAy1GSEhGRKstagxoCPFGmfCEwrHrhiIiIhKw1qLnAcUQfVLHjgT9lfTEz2x2YBuxGJLcj3H1uhe0HEUPbb3X3b2Z9HRGpTOtBSSPImqD+C7jLzPZj3b1Q+wLbA/+ZZQdmtiFwC3ABMXXSOGCWmW3n7p3N73c28E7g1oxxivQ71U4mWg9KGkWmJj53nwfsAfyKGCAxjEgaO7n7HzK+1n7AIHe/wN1XuftM4GHgsHIbp2T4IeCXGfcv0u8Ukkn7kqXQ0bE2mSyePafH+9R6UNIoMq+o6+6PAV/uxWuNAh4pKXsU2LV0QzPbDLgCOJioRYk0pUrJpKe1Ha0HJY2i0wRlZvcCH3b3l8xsLtDR2bYZ14MaApTO2bcCGFxm22nAJe7+kJll2LVI/5RHMtF6UNIoKtWgfgYULt1uq8JrvQJsVFI2GFheXGBmE4gmxAuq8JoiDS2PZKL1oKRRdJqg3P2Mol/vBP7g7quKtzGzVuJm3SzmA5NLynYCflBSdjiwN/BSqj29GfhPM9vL3T+W8bVE+oU8konWg5JGkbUP6k5gS2BJSfmOwHW8sWbU2T5a0pIdFxOj+HYDflK8kbt/uPh3M7sZuE/DzKUZ5ZVMtB6UNIJKfVDHAIVaVAsw38xK+6GGAH/J8kLuvtLMDiD6l84EngLGuvsSM2sDLnP3Id2MX6TfUzKRZlWpBnUF0W80ALgaOIv1V9XtIPqP7sj6Yu7+EPDeMuUzgLJjXN19bNb9i4hI/1GpD+p1Uv+QmT0J/B7Y2N1fTGXvBv6imcxFRCQPWefiWwI48LWisluBB81sZNWjEhGRppc1QV0E/JZ1fVIAOwD3EAMeREREqiprgtobONPd196z5O6vErM8/HsegYmISHPLmqBeBHYpU/52QMu+i4hI1WW9D+pK4HIz2xaYR4zg2wM4Dbgqp9hERKSJZU1QU9K2pwFbpLLFwPnAd3KIS0REmlymBOXua4BvAN8ws2HAygprOImIiPRa5uU2zOxdxJIZA9PvLUArsKe7T8onPBERaVaZEpSZnULMJLGcmLz1n8C/pIdvzyc0EQEtzy7NK+sovknAV9x9E+BvxCSvbyGWf5+bU2wiTS+PFXVFGkXWBLUlcFP6+T7gPe7+PPBVYHwegYmIlmeX5tadqY4KK6Q9Buyefv4rsHW1gxKRoOXZpZllTVC3EPdBvZNY1+mzZjYGOAl4Oq/gRJpdZyvnanl2aQZZE9SXib6mXYjl339DLLPx2fSYiORgxPg2BrS2rlem5dmlWWQdZj4OOMXdC+0KE8zseOA1Lbchkh8tzy7NLGuCuhC4F1jb8F08cayI5Ecr6kqzytrEdw9wcJ6BiIiIFMtag1oDnGNmpwJPAq8WP+jue1c7MBERaW5ZE9Q96Z+I1JhmkpBm1WmCMrMBaZJY3P2MzrYTkfwUZpIo3KxbmEkCUJKSfq9SH9QqMxteXGBmo82stbMniEh1aSYJaWaVElRLmbLbiDn4RKQGNJOENLOso/gKyiUtEcmJZpKQZtbdBCUiNaSZJKSZZV6wUERqTzNJSDPrKkFNMLPiGSM2AD5jZkuLN3L3S6oemYgAmklCmlelBLUIOKak7Hng8yVlHYASlIiIVFWnCcrdt69hHCIiIuvRIAkREalLSlAiIlKXlKBERKQu1XSYuZntDkwDdgMWAke4+9wy2+0JnJ+2exm4EjjL3TtqGK6IiPShmtWgzGxD4BbgBmBTYAowy8w2KdluMPAz4EfAUGB/YAJwVK1iFRGRvlfLJr79gEHufoG7r3L3mcDDwGEl220L/MHdL3b31e7+OHAz8N4axioiIn2slk18o4BHSsoeBXYtLnB3p2j13lTzOgC4PO8ARUSkftSyBjUEWFFStgIY3NkT0tIe16ftpuUXmoiI1Jta1qBeATYqKRsMLC+zLWa2JXATsdz8B9391XLbiYhI/1TLGtR8wErKdkrl6zGzUcBc4AkiOb2Uf3giIlJPalmDuhNoMbPJwMXAOGIY+U+KNzKzzYBZwEx3/3IN4xMRkTpSsxqUu68kBjuMA14ETgHGuvsSM2srmjV9PLFq7zFmtrzo3/W1ilVERPpeTW/UdfeHKDNc3N1nADPSzxcCF9YyLhERqT+a6khEROqSEpSIiNQlJSgREalLSlAiIlKXajpIQqQZLJ49h0XTZ9C+9AVahw1lxPg2ho8Z3ddhiTQcJSiRKlo8ew6PX3ARrFkDQPuSpfE7KEmJdJOa+ESqaMEll61NTmutWRPlItItSlAiVbTmtde6VS4inVOCEhGRuqQEJSIidUkJSkRE6pISlEgVtW4xrFvlItI5JSiRKhoxvo0Bra3rlQ1obWXE+LY+ikikcek+KJEqKtzrpBt1RXpPCUqkyoaPGa2EJFIFauITEZG6pAQlIiJ1SQlKRETqkhKUiIjUJSUoERGpSxrFJ01NazeJ1C8lKGlai2fPYcHUaaxpbwdi7aYFU6cBWrtJpB6oiU+a1qLpM9Ymp4I17e0smj6jjyISkWJKUNK02pe+0K1yEaktJShpWhsMGdKtchGpLSUoaVoddHSrXERqS4MkpGFUe8Td6uWvdKtcRGpLCUoaQh4j7jYYMoTXly0rWy4ifU9NfNIQ8hhxpyY+kfqmGpTkpppNcu1LlnarPAs18YnUNyUoAarfv7N49hyeuHAqHa+/DkQieeLCqUAPm+QGDIA1a8qX91DrsKFlE1zrsKE93qeIVI+a+BrQ4tlzmDdxEnePPZR5EyexePacXu/viQunxsm6o2NtMunNfp+84uq1yamg4/XXefKKq3u2w3LJqVJ5BpvutWe3ykWktmpagzKz3YFpwG7AQuAId59bZrsRwFXAvsBi4AR3v70WMeZRk6jrmgmVk0lP91lu8EGl8i7lUIP6x7w/datcRGqrZjUoM9sQuAW4AdgUmALMMrNNymw+E3gAGAocBcw0sx3zjrEwUqy4JrFg6rQe1ySqvT/IoWZCDskkDznUoDSThEh9q2UT337AIHe/wN1XuftM4GHgsOKNzOztwF7A6e6+0t1/A/wUODLvAKs9UiyPkWcNkUxy0LrFsG6VZ9pnJ31N6oMSqQ+1TFCjgEdKyh4Fdi2z3SJ3f6WL7aqu2lfUjXKFPnDj8vf9dFbeF0aMb2NAa+t6ZQNaWxkxvq2u9iki1VPLBDUEWFFStgIY3MPtqq7aV9R5XKHnkUx2POpIGDiwZIcDo7yHql3jGT5mNCOPOzqe39JC6xbDGHnc0b3qz8tjnyJSPbUcJPEKsFFJ2WBgeQ+3q7oR49vWm60AendFXe39QSSTx793Maxeva6wl8mkcEKu5mCOPN778DGjq5488tiniFRHLRPUfGBySdlOwA/KbDfCzDZy91eLtpufc3xVP1HnceLPY5+F/VbzRJ1XnCLSPFo6OmozrUsaxbcA+C5wMTAOuBwY6e5LSra9F/gd8DXg34hBEu9x94cyvM72wJN33HEH22yzTVXfg4iIVNezzz7L/vvvD7CDuz9V/FjN+qDcfSVwAJGYXgROAca6+xIzazOz4ia8ccA7iHugrgSOzJKcRESk/6jpjbopyby3TPkMYEbR788QyUxERJqUpjoSEZG6pAQlIiJ1qT/OZj4Q4Pnnn+/rOEREpAtF5+qBpY/1xwS1FUBbm2YDEBFpIFsRI73X6o8Jai7wPuBvwOouthURkb41kEhOb1jZomb3QYmIiHSHBkmIiEhdUoISEZG6pAQlIiJ1SQlKRETqkhKUiIjUJSUoERGpS0pQIiJSl5SgRESkLvXHmSTqgpl9EjiHuEN6NjDB3Rd3su1OwCXA3sRaWf/r7lPrLc6i5+wA3Acc5O531VOMZnY4cBqwDbAIOM3df5JTXLsD04DdgIXAEe7+hrvhzWwEcBWwL7HG2QnufnseMfUixj2B89N2LxPrsJ3l7jW5kz9rnEXbDwL+CNzq7t+spxjNbGPgIuBAoAO4ETje3VfVWZwGXArsASwDprn7lFrEmJVqUDkws1HECWkCMBR4HJjZybYbAb8A7gD+hTiov2Vm/1ZPcRY9ZyAwHRiSd3zp9brzWb6b+GIeRXyWXwF+aGY75xDXhsAtwA3ApsAUYJaZbVJm85nAAyn+o4CZZrZjtWPqaYxmNhj4GfCjFOP+xOd9VN4xdifOEmcD76xBeEC3Y7w6bbM9sfDqXsSxWG9xzgB+DWwOfAD4opkdWIs4s1KCysdniCu737n7a8TS9f9uZm8rs+2BwD/dfYq7r3b3+4ia1KN1FmfBqcCDxBVXLXQnxu2A77n73e7e4e4/BxzYJ4e49gMGufsF7r7K3WcCDwOHFW9kZm8nTlCnu/tKd/8N8FPgyBxi6lGMwLbAH9z94nQMPg7cTJnFRfs4TgDMbD/gQ8AvaxQfZP97bwUcBBzl7i+nmv5BFC3IWg9xFsJN/7cQNb0O4LWaRJmRmvh6KF2pbF7moQ5gFDCvUODuK8zsGWBXogZQbE/gYTObBhwCvACcnVYZrqc4MbN9gU8C7wYOr0Z81YzR3W8kmlMK+31rev591Yq1yCjgkZKyR1NcpdstcvdXSrbbO4eYSmWK0d0dOLjwe/p7HABcnneASdbPEjPbDLiCiPfs/ENbK2uM7yKaltvM7IvAIOCHRLNzLWT+LIGziBrWmcSEree7+6x8w+se1aB67t+IGdNL//2VaP5aUbL9CmBwmf1sDnyCOAm/BTgGuMzMqnX1WpU4zWwIcC3weXcvfU5dxFgS7zbArcDV7v7nKsdLN+LqUfxV0u3XNrNW4Pq03bT8QltPd+KcBlzi7g/lHtX6ssa4OdG0twvRBzSGaCX5as7xFXTns+wATk7PeSdwiJnVomafmWpQPZQGB7SUe8zMbgE2KikeDCwvs3k78KC7X5l+v8vMfgKMBX5XR3FeDFzn7vf2NqZSVYyx8Jx9gZuAHwMnVifKN3glY1xZt8tDt17bzLYkPrc1wAfd/dV8w1srU5xmNgEYBlxQm7DWk/WzbCdqIye7+3JguZl9l7jwPCf3KLN/lnsBk91961R0v5l9GziW6POtC6pB5WM+69p3C53QI1J5qUeBzUrKanXh0J04DwNONrN/mNk/iEEIt5nZf9dRjJjZJ4iO32+5+wnuvqYWcSU7lYlrPjAiDYaptF0essZYGIwyF3iCSE4v5R/eWlnjPJxoGn0pHYMfBf7bzG7LP8TMMRb6jjctKqtlRSBrnNsCG5pZ8YXh60BNRhpmpRpUPq4Dfpc6c/8AfAv4i7s/VmbbG4EpZvZfwLlEx/SBRGdn3cTp7utdlaUTxNgaDDPPHKOZvYcYYfgpd78557juBFrMbDJRuxxHNOmsN6Td3d3M7if+xl8jmjMPAt6Tc3yZY0z9OrOAme7+5RrEVSrrZ/nh4t/N7GbgvhoNM88a44NmNg8438zGEzW+yUS/WS1kihO4m6jpnWFmZwI7AF8GanJ7S1aqQeXA3R8EjiDay5cCOxP9TACY2TQz+3na9m9EO/VHiAESVwNHV7oHpC/i7CvdjPErwIbE0PLlRf+OzSGulcRAgnHEvWunEAl7iZm1mVlxk8o4YrjxYuL+oiNr0YfSjRjHk/o/Sz636/OOsZtx9pluxvgRYjTcQqJv+RZq1CyZNc40uvAA4P3E92oWcA1x/1bd0Iq6IiJSl1SDEhGRuqQEJSIidUkJSkRE6pISlIiI1CUlKBERqUtKUCIiUpd0o24FZvYUMUN2OVPc/dTaRdNz6X2c6+4XZ9j2zcRU/R8E5rr7+3rxui3ARGC6u7+Wpqo5192H9XSf1VQaXy/2MxY42N0/Z2bXAkPc/dCMzz0HOJ6YF2373szgYGZjiJnx7zOz7YEngV37YN66hmdmHcDH3b0Ws1TkKs2juQx4f29vrDezQ4GPuvvnqxFbV1SD6trXiYXySv/9T18G1U3vJm4AzmIskZz+nZi5vDdGEzNiFy6EbiBmW64XpfF1W1qc7rvA6anoRCLpZXnuW4jlQ04Gdq/C9EJ3EQs1AjxDHKe1WLalP9oK+FVfB1Fv0ooBu6SZXXKnGlTXlrn7830dRG+4+5JubL4p8Hd3/+oHZQkAAA3KSURBVFMVXnq9CWDT5KO1moA0i7IT1HbTF4A/u/vTAO7+z248tzBf26/d/akqxLKWu68GGvq47UuN/p3P2UXAN6nBdGxKUL2Qlib4M/CEux+Uyk4mFvXb1d2fNbP3E7MY704sIfE/7n5F2vbtwIXElfwSooZxmru3p8dPI06AWxBrvHw9LcKHmX2BmMJ/W6Ip5xx3/0EncT5FauJLTVAriBmPDyWmV7rS3c82s28C30jP6SCW1rjWzD6b3tM2xBX5ae7+s6L9H0PMN7YN8BBwEvAsMS8YwDIzKzQJrG3iS4sOnpve/xpi9vGT3H1ZURPVJ4h1a7YlJjM9Oq1fVPoe9yPmNbwmfWYz3X1SmpPsGGIJhOXEyrHHEHOkrRdfeq8HEPP9GTFVzbnufk25zzU5nphOphDHtaQmvtSkeTyxUu1kYm2gXwCTiFpt4fUXmtn33X2Cme1D1Mj2JGpBV6QY1qT9lz2e0t8Y4FYz+z5xAlnbxGdm/0Ksn3QIMTnxHODEwmeZnn8+sc7S3sSksae4+63p8YPT3+Gt6XUvcffvdPahpOPnbOJv8TyxrPhIKh/vuwPnAfsC/yCWID87PbYV8G3gw8Cb0ud4orv/zcymA29290OKXv9LwLHu/naL5eHPIlYJ3ohYKv6LJe/9R8CniIuW3Yhpgj7u7rdleP77Uty7Ai8R6z99LV0klH4u1xIzng9Jn/Vz6TO4vmibTr9v6fkbAG9L/z7p7r8ueY3B6XP+BNG0d2rJ4/9KHGP/QUz6/AwxufKVZvZJYlmdf3X3ZWn7EcSx9DZ3X0gsY3Olme3u7veXvsdqUhNfL6Qv1hHAx8zsIDMz4kt5QkpOOxFfpN8S662cClxsZh8yszcRK4IuJBY5Gw/8J/A9WHtC+BKxouxOxIn1/8xsEzPbg7iKOQkoJLlrrfJKuMUmEglkL2JuuLPMbE8iWXw9PbYVcIOZfTjt/3TiC3gZcGOamJWUeM4jmjx3I066PyO+GOPS640kTkZrmdnmxHIiq4D3ESfO9/LGpshvEif096eYzq3wvoamz2MP4FwzOzw9/yTiyzyBmKh1EvGlXC8+i6XhbyLm/duFWMjtPDP7VLkXS9tvT+WVXXdL7++DxBLqhwBHA79n3RXo3sCJZjY87esXxGf9ReA40lpClY4nIuFBHEfllhi5kfgMDydWGH6NWAq8eJ2gM9N734tYDPIaM9swndBuIBKYEXMeTjGz/Su8b4A2YinxzxHJudLxPgz4DXHC3ps4Rk82s4kpQdxBzGL/kbTPtwA3p37E64ADUnNrwadYt4rtmcTM54el9+7A7JS0C44iLtgOLtPU2unzzWwgsfrwbcR8i59N+5pQ4XOZALxMHKcXE3NH7pc+h4rft6LP9dL0Ofy+zP4vIY65jxDH+OSSx38ADAf2J5rcfwpcarHcyk+J7+RBRdt/GrgnJSfS53MvMZdfrlSD6tp5Zlauv2mUuy9y93vM7AJiMsjngNvd/YdpmyOBh929sFjZYxYzR0OcKFYBx7l7B+BmdjTwWzP7KnHiaweedvenLGYcnpOesx1R43g6NS1dYmaPE1elWSxw98IKn2enWsae7v4nM1sGrC40cZjZ14HveCwdDbAgJbOTiS/0scBl7n512r6w/MZmxFUowGJ3fzXy91qfJi6QxhfWHUo1jj+mmuXKQnzuPic9fgnQ1fIe33L3BWn7rYEJRR3dT5vZbOJvt9rMSuP7KjDD3QsL9S0ws5HELM8zeaO9gOfd/YUK8Qwilv9+nlg5+RfEZ73SzArPW+Lu/zSzM4iBKWel8sctZkC/kLgA6PR4SpOBAvwj7WvtEi5mtguRIN/t7vNSWRvwNHGyK8y0fX3h75xiuZ84Dt+c3scz6Xh72sz+DpSbnb/Y5e7+cNrf56l8vB+WHp+YJjydbzHJ72qi1vRWYhmQ59L+DiOu6j9I9BW9DHwcuM7MtiOS3HiLZU4mA/u7+90pri+mRDCeSBAAN3iZtc4yPP86YpHCv6Vm2qfM7D+IiYE783T6HNYAj6bkdAzRh9jV9w1ikvxry+3YzDYhvlsHF+I1s0nESgAFPwNuKyQcMzubuBh+m7v/1sxuIhJ84Tz2aSJRFptPHP+5UoLq2rdY94cq9lzRz6cSV8a7EoMMCtZbrhzA3S8BMLNzgR2J5qXCwy3ESftt6TU/RzT//Im4srkmnUh/QUyXf5+ZPUxcvV3j7v/I+J5Kl3NfRpyAytkZ2CedKAsGse7kNIq4si68vw7gv9J73IbO7UwslVDcJzWXSEzFy7QXnwRfrhBnwYKiWGab2V7pC7hTek0jriA7i2nXVPMq2IDO18j5V2Im6EpK+zBfJk74nb3+frb+zNgDgI3MbCgVjqcu7Ex8rmv7Fd39FTP7S3qsoPSzhvi87yNqI780s4XECe4H7v731LxVPOP9D9396PTzgqLynal8vI8CHkjJqRDj9QAWS9E8XUhO6bFnU9Pczu7+KzP7P2JQz3XEyXWeuz+eknMrUVssnhn7Tay/blJxrMVGVnq+u79oZt8CLrdoIr+dSPSVViP4va+/Ttm9xIUCdP19qxQrxHsaRHQ9FMwjLmgLLgU+YWaFFph3pfKB6f/pxN96c2Bromb4o5LXeYH4TuVKCaprS939iS622Zo4WQ0irtwK/TMrieHD5WxAXNWUG67515SI9iSq4R8lmgUmm9no1KfwIWJtoY8DHyOaiD7m7ndkeE8ry5R1NmBgA2Kk2a0l5YWTdqX3WElngyUKJ62C0li7Gtiwdr+pRnYp0S/1c6If4YwKz92AaDrNutT5mgzxdPezvomSPoPkn9T+swZoSRcdn0kXVQcRTUfHplrRjURzY8HLRT8Xv27F453K7y1L/DOAO1Mz32GsuwgpnOM+xBtrNZ3FWqzL57v711PfUOGz+ZWZnebuna2g+3rJ7wOImmLh9Sp93yrFWqz4GFtd2H9qEv0lcc6aSTSdPkw0WxbcBfyd6CPbEZjlbxxoNbAo5twoQfVS+oNfRfzRHwQuM7OdPUZzPUa0ExdvfwVxYD9CXOk96+keHDPbm+hvONKis37LdIU8y2LwxQLgI6m9/r2pA/luYlXRu4n25iwJqjseAbYrTtKp2W8N0ez0GNGWXtzJ+zDRhl6pdvEIcdLbqKgWtReR5Ks1NPo4ornk9BRXC3G1XqiFlJ4QHwHeWvJeJxJX9yeV2f/zxACWankEOKDk9T9OdHZPoMLx5O4nd7HfDYmBF4UmvsHEQIsfdxWUxeCFz7n7SURt6gwzuw443N2nEwMqsry3To/39N4OM7NB7r4qPX5aivEyYDsz28pj/bRC8+12pGPF3f9gZs8R/Yu7sa5J9gkiIQx399+n57YQtYSZROtDJRWfb2YPEAnlS2nQyHcs7m37DJ0v8b5Hye97E82phc+p0vetK48SyX4f1i1SuCvrWh5GEX2RI4ua+PZOj7VAtIKkv+/HiQRV7nWHUYNRokpQXds4dR6WWunuLxId3nsR1eAlRHvtuURH6SVEzeYs4PtER/Z44irrHuIk/v30+MbEgIVFqQ9hIHGw/504qewLbJl+XgF8Iz02i6hqvyM9v9q+DVxvZo8Sbf0fIGoh49Pj5xMjeu4jRjdNJGqTc4hObYA9UzNlsRnAacD01DSyGfF5/drd51uM4uutF4D3Wyxn3kL0JYwiLiQgRvUVx3cucE86IfyIODmeTwx8KedPwHAz29KrMyx5KtG/cWH6eVvi5HyLu69JfXCdHU+F97OLma3XcZ6aun5MDHo4lhhp9g3iCrhc31qpF4nFDAsj1LYmjsfp3XhvP6Ty8T6DOK4uKWr+nkxcGPyaOIHPTM1SEKPQHkuPFVyX3tedhb+Huy9Pn9v3zKydaN6eTDTFf6OroDM8fwmpb8jMzgM2IWpblZr43pWanX9AtH4cQNx3CF1/37qKd5mZXQl8N/29Xia+V4WLsX8Qf/dPpc+8MMgKoimzYHp6D6uJBRdL7c66QSi50Si+rp1DDKst/fdjM9sW+F/gG+7+TLoyPI6oAX0wdSh/nGiie5gYUTbR3X/j7q8Qwzw3I9qgbyUOiE8DuPv/EV+A7xBfxLOB49Nz7yWuOk8iquZXAOd55eHQPeLuPwFOSK81nxjBdWyhEzf9fyowhTjxjyZqAUvS77cRSfQLJftdQYzi2iS97x8To/oOoXpOJL6Y84gTWSvRp1i4gl0vPo97vw4lmogeJk6C/0OcNN7A3R8hOurHVCNYd3+WGBCwF3FCnk6MnvtSerzT4ynt4lwiCVxVZvdHEMfZT4kLicHA6HSR1VVczxB/l4PT695EjFyb0o331tXx/jJxot4pvfdpxGwt16YmxrFEMriLaCV4jhi4UNwsOYMYvn1dyct/lXW3IDxA/P0PKAymyaDT56fa/8eIWsp9xHH2EDECszOziJr8/cTf5eDCAI2uvm8ZTSb+zj8mWna+T2q+dfe/Et/FLxC1rQuJBPYAUcMmbfcQkYxvTn+7tSxGP+5G17XPXtOKuiK9YGZfAd7n7gf2dSxS/6ybU2H1ldSC8ywxCvaXJY9NIu6/6uo2g15TDUqkdy4F9jCzHfs6EJFqMLNxRBPzckqme0r9b0fTebN3VSlBifSCuy8nmuDO7OtYRKrkLKLZ8vMlw+EhmsAfdPc73/i06lMTn4iI1CXVoEREpC4pQYmISF1SghIRkbqkBCUiInVJCUpEROrS/wOhdE6oYpEMmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "plot_sweep_frame_difference(frame)\n",
    "\n",
    "decorate(xlabel='Excess infection rate (infections-recoveries per day)',\n",
    "         ylabel='Fraction infected',\n",
    "         legend=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# The results don't fall on a line, which means that if we know the difference between\n",
    "# `beta` and `gamma`, but not their ratio, that's not enough to predict the fraction infected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose you run a survey at the end of the semester and find that 26% of students had the Freshman Plague at some point.\n",
    "\n",
    "What is your best estimate of `c`?\n",
    "\n",
    "Hint: if you print `frac_infected_series`, you can read off the answer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.211261    0.999900\n",
       "4.642296    0.989902\n",
       "3.987365    0.979904\n",
       "3.612133    0.969906\n",
       "3.350924    0.959908\n",
       "3.151808    0.949910\n",
       "2.991711    0.939912\n",
       "2.858363    0.929914\n",
       "2.744467    0.919916\n",
       "2.645332    0.909918\n",
       "2.557767    0.899920\n",
       "2.479505    0.889922\n",
       "2.408879    0.879924\n",
       "2.344627    0.869926\n",
       "2.285771    0.859928\n",
       "2.231541    0.849930\n",
       "2.181315    0.839932\n",
       "2.134590    0.829934\n",
       "2.090947    0.819936\n",
       "2.050040    0.809938\n",
       "2.011573    0.799940\n",
       "1.975299    0.789942\n",
       "1.941002    0.779944\n",
       "1.908499    0.769946\n",
       "1.877628    0.759948\n",
       "1.848249    0.749950\n",
       "1.820238    0.739952\n",
       "1.793487    0.729954\n",
       "1.767898    0.719956\n",
       "1.743384    0.709958\n",
       "              ...   \n",
       "1.181034    0.290042\n",
       "1.173263    0.280044\n",
       "1.165630    0.270046\n",
       "1.158132    0.260048\n",
       "1.150765    0.250050\n",
       "1.143524    0.240052\n",
       "1.136407    0.230054\n",
       "1.129409    0.220056\n",
       "1.122527    0.210058\n",
       "1.115758    0.200060\n",
       "1.109099    0.190062\n",
       "1.102547    0.180064\n",
       "1.096099    0.170066\n",
       "1.089751    0.160068\n",
       "1.083503    0.150070\n",
       "1.077350    0.140072\n",
       "1.071291    0.130074\n",
       "1.065323    0.120076\n",
       "1.059444    0.110078\n",
       "1.053651    0.100080\n",
       "1.047943    0.090082\n",
       "1.042317    0.080084\n",
       "1.036772    0.070086\n",
       "1.031305    0.060088\n",
       "1.025914    0.050090\n",
       "1.020598    0.040092\n",
       "1.015356    0.030094\n",
       "1.010185    0.020096\n",
       "1.005083    0.010098\n",
       "1.000050    0.000100\n",
       "Length: 101, dtype: float64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "frac_infected_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# It looks like the fraction infected is 0.26 when the contact number is about 1.16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
